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Introduction 


1.1  Background 

Multi-anchored  or  tieback  wall  systems  are  often  used  for  temporary  support 
of  excavations  that  have  space  restrictions  due  to  adjacent  structures,  highways, 
railroads,  etc.  In  some  cases,  multi-anchored  systems  may  remain  as  permanent 
structures  after  construction.  In  U.S.  Army  Corps  of  Engineers  projects,  perma¬ 
nent  tieback  wall  systems  are  used  as  guide  walls  and  approach  walls  on  naviga¬ 
tion  projects,  and  as  retaining  walls  on  highway  and  railroad  protection  and 
relocation  projects. 

The  behavior  of  multi-anchored  systems  may  be  strongly  influenced  by 
factors  such  as  the  sequence  of  excavation  and  installation  of  anchors,  fluctu¬ 
ations  in  the  water  table,  and  the  nonlinear  stress-strain  behavior  of  soils. 
Therefore,  to  obtain  accurate  predictions  of  the  magnitudes  of  stresses  and 
deformations  in  the  structure  and  the  surrounding  soil,  it  is  necessary  to  perform 
soil-structure  interaction  (SSI)  analyses  that  model  the  construction  and  operation 
stages  of  the  system.  For  such  analyses,  adequate  modeling  of  the  excavation 
process  is  required. 

A  substantial  amount  of  research  has  been  performed  in  recent  years  on 
another  type  of  earth-retaining  structure:  lock  walls  for  navigation.  These  studies 
have  included  SSI  analyses  of  the  Red  River  Lock  and  Dam  No.l  (Ebeling  et  al. 
1993;  Ebeling  and  Mosher  1996;  Ebeling,  Peters,  and  Mosher  1997),  the  North 
Lock  Wall  at  McAlpine  Locks  (Ebeling  and  Wahl  1997),  and  Locks  27  (Ebeling, 
Pace,  and  Morrison  1997).  These  are  good  examples  of  state-of-the-art  tech¬ 
niques  available  for  SSI  analyses. 


1.1.1  Common  types  of  multi-anchored  systems 

Multi-anchored  systems  can  be  constructed  using  different  materials  and 
configurations.  The  following  are  the  most  common  types  found  in  practice: 

•  Vertical  sheet-pile  systems  with  wales  and  post-tensioned  tieback 
anchors. 

•  Soldier  beam  systems  with  wood  or  reinforced  concrete  lagging  and 
post-tensioned  tieback  anchors. 
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•  Secant  cylinder  pile  systems  with  post-tensioned  tieback  anchors. 

•  Continuous  reinforced  concrete  slurry  wall  systems  with  post-tensioned 
tieback  anchors. 

•  Discrete  concrete  slurry  wall  systems  (soldier  beams  with  concrete 
lagging)  with  post-tensioned  tieback  anchors. 

Figure  1.1  illustrates  the  use  of  a  multi-anchored  system  for  a  typical  navi¬ 
gation  project.  Because  of  the  space  restrictions  imposed  by  an  adjacent  railroad, 
excavation  for  the  expansion  of  the  waterway  requires  use  of  a  multi-anchored 
system.  For  simplicity,  it  is  assumed  that  the  multi-anchored  system  depicted  in 
the  figure  corresponds  to  a  continuous,  reinforced  concrete  slurry  wall  with 
tieback  anchors.  Tiebacks  consist  of  post-tensioned  tendons  with  a  grouted 
anchor  region.  A  berm  of  granular  material  or  riprap  is  placed  at  the  toe  of  the 
wall  to  minimize  erosion  and  improve  stability. 


Figure  1 .2  illustrates  the  typical  construction  sequence  of  a  reinforced  con¬ 
crete  slurry  wall.  Initially,  a  trench  is  excavated  using  a  clamshell-type  tool.  The 
excavation  is  stabilized  by  the  use  of  mud  slurry.  The  finished  trench  acts  as 
formwork  for  the  reinforced  concrete  panel.  Placement  of  the  concrete  using  a 
tremie  pipe  displaces  the  mud  slurry  and  leaves  a  structural  concrete  wall  that  can 
be  excavated  and  tied  back  in  much  the  same  manner  as  the  other  tieback  wall 
systems.  The  walls  are  reinforced  using  preassembled  cages,  which  are  dropped 
into  the  slurry  trench  just  before  concrete  placement.  Slurry  wall  systems  are 
usually  0.6  to  0.9  m  (2  to  3  ft)  thick  and  can  be  placed  to  depths  of  30  m  (100  ft) 
or  more.  The  construction  process  can  be  summarized  as  follows: 
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CONCRETE  BLURRY  LEVEL- 

PANEL—  1 


GROUND  LEVEL 


CONCRETE 
PANEL- 


STOP  END 
TUBE 


GROUND  LEVEL 


a.  Excavation  under  mud  slurry 


b.  Placement  of  stop  end  tube 


CONCRETE 
PANEL- 


REINFORCEMENT 

CAGE 


.STOP  END 
TUBE 


TREMIE  PIPE 


—.STOP  END 
1/  TUBE 


GROUND  LEVEL 


c.  Placement  of  reinforcement  cage 


d.  Placement  of  concrete 


_  Figure  1 .2.  Typical  construction  sequence  of  a  reinforced  concrete  slurry  wall  (from  Strom  and  Ebeling 

•  2001) 

a.  Guide  walls  are  constructed  to  facilitate  positioning  and  alignment  of  the 
clamshell  during  the  excavation  process.  To  stabilize  the  excavation, 
mud  slurry  is  kept  inside  the  excavation  to  a  level  above  the  water  table. 

•  As  illustrated  in  Figure  1 .2a,  the  excavation  for  each  panel  follows  a 

staggered  sequence.  Two  end  excavations  are  performed  first,  leaving  a 
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central  core  intact.  After  the  end  excavations  are  completed,  the  central 
core  is  removed. 

b.  A  stop  end  tube  is  placed  at  one  end  of  the  panel  excavation.  This  tube  is 
extracted  after  concrete  placement,  leaving  a  semicircular  indentation. 
This  indentation  serves  as  a  guide  for  the  excavation  of  the  adjacent 
panel  and  allows  the  creation  of  a  shear  key  between  the  panels. 

c.  Once  the  panel  has  been  excavated  to  the  desired  depth  and  the  slurry 
cleaned  of  fine  excavation  material  (desanded),  the  reinforcement  cage  is 
lowered  into  the  excavation. 

d.  One  or  more  tremie  pipes  are  used  to  place  the  concrete  without 
contamination  from  the  slurry. 

e.  Once  the  wall  is  finished  and  the  concrete  reaches  its  desired  strength, 
the  excavation  and  tieback  installation  process  can  begin. 

Construction  of  the  second  navigation  lock  at  Bonneville  Lock  and  Dam 
required  the  use  of  concrete  slurry  walls  to  retain  the  foundation  of  an  adjacent 
railroad  line.  Detailed  descriptions  of  construction  procedures  for  the  continuous 
reinforced  concrete  sluriy  wall  and  for  the  discrete  sluny  wall  systems  used  at 
Bonneville  Lock  are  presented  by  Munger,  Jones,  and  Johnson  (1990)  and 
Maurseth  and  Sedey  ( 1 992),  respectively. 


1.1.2  Response  of  multi-anchored  systems  to  excavation 

A  waterway  expansion  project,  such  as  that  presented  in  Figure  1.1,  requires 
SSI  analyses  to  determine  the  magnitude  of  the  deformations  of  the  soil  above 
the  excavation,  and  the  bending  moments  and  stresses  in  the  retaining  wall.  Such 
analyses  require  close  modeling  of  the  construction  stages  of  the  multi-anchored 
system.  The  finite  element  analyses  performed  by  Mosher  and  Knowles  (1990) 
for  the  tieback  walls  at  Bonneville  Lock  and  Dam  are  a  good  example  of  the 
available  techniques  that  can  be  used  in  SSI  analyses  of  multi-anchored  systems. 

Figure  1.3  illustrates  some  of  the  construction  and  operation  stages  of  the 
hypothetical  navigation  project  shown  in  Figure  1.1.  For  simplicity,  it  is  assumed 
that  construction  is  performed  in  the  diy.  After  completion  of  the  continuous 
reinforced  concrete  sluny  wall  (Figure  1.3b),  the  soil  in  front  of  the  wall  is  exca¬ 
vated  to  an  elevation  slightly  below  the  position  of  the  first  row  of  anchors.  The 
anchors  are  then  installed  and  tensioned  according  to  the  project  specifications 
(Figure  1 .3c).  Once  these  anchors  are  tensioned  and  tested,  excavation  continues 
until  reaching  the  position  of  the  second  row  of  anchors.  The  process  is  repeated 
until  reaching  the  bottom  of  the  excavation  (Figures  1 .3d  and  1 .3e).  Once  the 
excavation  is  completed,  the  granular  toe  berm  is  placed  against  the  toe  of  the 
wall  (Figure  1 .3f).  During  operation  of  the  navigation  facility,  the  water  level 
outside  the  wall  reaches  its  normal  elevation,  which  may  fluctuate  periodically 
during  the  life  of  the  structure  (Figure  1.3g). 
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a.  Initial  state 


e.  Final  excavation  and  installation  of 
bottom  row  of  anchors 


b.  Construction  of  reinforced  concrete 
slurry  wall 


f.  Construction  of  riprap  blanket 


c.  First  stage  of  excavation  and 
installation  of  anchors 


g.  Operational  stage 


d.  Intermediate  stage  of  excavation  and 
installation  of  anchors 
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1.1.3  Need  for  nonlinear  finite  element  SSI  analyses 

Excavation  in  soils  is  probably  one  of  the  few  geotechnical  problems  that 
cannot  be  solved  analytically,  except  in  tunneling  where  the  excavated  region  can 
be  approximated  by  a  circle.  There  are  also  analytical  methods  to  determine  the 
stability  of  natural  slopes  and  estimate  the  magnitude  of  lateral  earth  pressures  on 
structures.  However,  there  are  no  simple  and  closed-form  analytical  solutions  that 
can  be  used  to  estimate  or  predict  the  magnitudes  of  deformations  expected  from 
an  excavation  in  soil,  even  if  the  soil  is  assumed  to  behave  as  a  linear  elastic 
material.  The  situation  becomes  more  complicated  when  soil  nonlinearity,  the 
presence  of  water,  interaction  of  the  soils  with  structures,  and  construction  effects 
are  taken  into  consideration.  In  many  cases,  a  reliable  analysis  of  the  response  of 
soils  and  structures  due  to  construction  and  operational  loads  can  be  obtained 
only  by  performing  nonlinear  finite  element  (NLFEM)  soil-structure  interaction 
analyses. 

Two  methods  that  have  been  used  in  Corps  of  Engineers  projects  are  the 
beam  on  rigid  supports  (Strom  and  Ebeling  2001)  and  the  WINKLER  method. 
The  beam  on  rigid  support  method,  referred  to  as  the  RIGID  analysis  in  this 
report,  does  not  take  soil  deformations  into  account  and,  although  widely  used, 
has  limited  application  for  those  circumstances  in  which  actual  load- 
displacement  characteristics  of  the  system  are  desired.  The  WINKLER  method 
relates  soil  pressures  to  wall  deformations  and,  for  that  reason,  can  be  considered 
an  improved  method  of  analysis,  especially  when  used  in  a  staged  construction 
analysis.  However,  it  has  not  been  a  reliable  tool  for  evaluating  soil  displace¬ 
ments  that  occur  either  in  the  soil  mass  or  at  the  ground  surface  of  the  soil 
retained  by  the  tieback  wall.  Wall  displacement  information  provided  by  the 
WINKLER  analysis  is  unreliable  because  the  value  of  the  coefficient  of  hori¬ 
zontal  subgrade  reaction  is  dependent  on  the  extent  of  the  zone  of  influence,  a 
quantity  that  is  difficult  to  properly  establish.  In  addition,  arching  that  takes  place 
in  the  retained  soil  influences  the  soil  displacement-pressure  response,  but  this  is 
not  included  in  the  WINKLER  analysis. 

The  NLFEM  can  overcome  the  shortcomings  of  other  analysis  procedures, 
such  as  the  RIGID  and  the  WRINKLER  methods,  since  it  permits  a  complete  and 
precise  solution  based  on  the  stress-strain  laws  for  the  soils  involved,  the  bound¬ 
ary  conditions  of  the  problem,  and  the  basic  equations  of  mechanics.  Although 
the  finite  element  method  is  the  most  suitable  one,  the  complexities  of  the 
method  are  such  that  it  is  not  routinely  used  in  the  design  or  evaluation  of  tieback 
wall  systems. 

In  several  cases,  the  NLFEM  has  been  used  with  success  to  evaluate  the  per¬ 
formance  of  tieback  wall  systems.  In  particular,  application  with  respect  to  the 
performance  evaluation  of  a  temporary  tieback  wall  used  to  facilitate  construc¬ 
tion  of  the  Bonneville  navigation  lock  (Mosher  and  Knowles  1990)  is  described 
below.  The  NLFEM  provided  a  detailed  and  accurate  representation  of  the  tie- 
back  wall-soil  system  response  to  various  loadings  that  occurred  prior  to,  during, 
and  after  wall  construction.  The  objectives  of  the  NLFEM  were  to  confirm  by 
incremental  analysis  that  ground  movements  behind  the  wall  during  and  after 
construction  would  meet  stringent  displacement  requirements.  The  results  of  the 
analysis  were  confirmed  by  instrumentation.  The  Bonneville  navigation  lock 
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NLFEM  analysis  is  described  in  detail  to  demonstrate  the  differences  between 
actual  behavior  and  that  usually  assumed  in  the  design  of  tieback  wall  systems. 
These  differences  are  especially  noticeable  with  respect  to  stiff  wall  systems  with 
high  anchor  prestress  loads. 


1 .2  Review  of  Previous  Work  on  SSI  Analysis 

The  first  systematic  SSI  analyses  of  retaining  wall  behavior  were  presented 
by  Clough  and  Duncan  (1969,  1971)  and  Duncan  and  Clough  (1971).  These 
investigators  used  the  hyperbolic  constitutive  relationship  developed  by  Duncan 
and  Chang  (1970)  to  model  the  behavior  of  the  backfill,  and  extended  it  to  model 
the  behavior  of  the  wall-to-soil  interfaces.  Relative  movement  at  the  interfaces 
was  achieved  using  the  joint  element  developed  by  Goodman,  Taylor,  and 
Brekke  (1968). 

In  their  analyses  of  Port  Allen  and  Old  River  U-frame  locks,  Clough  and 
Duncan  (1969)  and  Duncan  and  Clough  (1971)  demonstrated  the  importance  of 
close  modeling  of  the  construction  stages  of  the  lock  and  backfill  placement. 
They  demonstrated  that  a  simple  linear  elastic  model  for  the  soil  and  use  of 
gravity  turn-on  analyses  are  not  adequate  to  model  the  behavior  of  the  soil-lock 
system.  They  also  proved  that  the  downdrag  or  vertical  shear  force  exerted  by  the 
backfill  on  the  wall  has  an  important  influence  on  the  behavior  of  U-frame  locks. 
Their  work  provided  fundamental  understanding  of  previously  unknown  aspects 
of  lock  wall  behavior. 

Clough  and  Duncan  (1971)  presented  a  systematic  approach  to  SSI  analyses 
of  retaining  wall  behavior.  They  observed  the  importance  of  modeling  the 
different  stages  of  construction  of  the  wall  and  placement  of  the  backfill  in  the 
SSI  analysis.  They  found  that  when  the  stages  of  placement  of  the  backfill  were 
closely  modeled,  the  resulting  horizontal  and  vertical  loads  acting  on  the  wall 
were  substantially  larger  than  those  obtained  using  classical  earth  pressure 
theories.  The  results  of  these  analyses  were  consistent  with  some  previous 
experimental  work  and  field  observations. 

Ebeling,  Duncan,  and  Clough  (1990)  performed  a  comparison  between 
results  from  conventional  equilibrium  and  finite  element  analyses  of  several 
hypothetical  gravity  walls  founded  on  rock.  Their  analyses  were  performed  with 
the  backfill  placement  analysis  option  incorporated  in  SOILSTRUCT  (Clough 
and  Duncan  1969).  A  range  of  possible  values  of  shear  stiffness  was  assumed  at 
the  interfaces  between  the  wall  and  the  backfill,  and  between  the  backfill  and  the 
rock.  Ebeling,  Duncan,  and  Clough  (1990)  concluded  that  the  magnitude  of  the 
downdrag  force  is  significantly  affected  by  the  concrete-to-backfill  and  rock-to- 
backfill  shear  stiffness  values.  They  also  concluded  that  the  conventional 
equilibrium  analyses  neglect  the  true  process  of  soil-structure  interaction  and 
tend  to  yield  very  conservative  results. 

Ebeling  et  al.  (1992)  performed  analysis  of  several  hypothetical  gravity  walls 
founded  on  rock.  The  hypothetical  walls  were  based  on  several  representative 
examples  of  lock  walls.  Ebeling  et  al.  (1992)  found  that  conventional  equilibrium 
analyses  are  very  conservative  because  they  do  not  account  for  the  stabilizing 
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effect  of  the  downdrag  forces  generated  by  the  settlement  of  the  backfill.  At  the 
time  of  their  work,  it  was  not  known  whether  these  vertical  shear  forces  persisted 
under  field  conditions  nor  whether  they  could  be  relied  upon  for  the  stability  of 
the  structure.  These  researchers  also  indicated  that  the  behavior  of  retaining 
structures  founded  on  soil  might  differ  substantially  from  that  of  structures 
founded  on  rock.  In  soil-founded  structures,  the  concrete-to-foundation  interface 
is  not  bonded  as  in  the  case  of  concrete-to-rock  interfaces,  and  relative  interface 
displacements  may  occur,  including  more  redistribution  of  the  earth  pressures. 

Experimental  data  from  the  Instrumented  Retaining  Wall  Facility  at  Virginia 
Tech  (Filz  1992)  showed  that  downdrag  forces  on  the  nonmoving  test  wall  were 
significant  and  tended  to  either  remain  constant  or  increase  with  time  after 
backfill  placement. 

Ebeling  et  al.  (1993),  Ebeling  and  Mosher  (1996),  and  Ebeling,  Peters,  and 
Mosher  (1997)  presented  the  results  of  extensive  SSI  analysis  for  the  soil- 
founded  Red  River  Lock  and  Dam  No.  1.  A  reinforced  soil  berm  was  recom¬ 
mended,  among  other  alternatives,  as  a  solution  to  problems  induced  by  siltation 
of  the  lock.  The  SSI  analysis  procedures  were  validated  against  instrumentation 
measurements  from  the  lock  taken  at  the  end  of  construction  and  several  opera¬ 
tional  stages.  Their  analysis  revealed  that  important  changes  in  normal  stresses 
may  occur  at  the  soil-to-structure  interface  during  backfill  placement  and 
operation  of  the  lock,  and  underscored  the  importance  of  selecting  appropriate 
interface  stiffness  values  for  these  loading  conditions.  They  also  noted  that 
conventional  equilibrium  analyses  are  inadequate  for  the  design  of  this  type  of 
structure. 

Ebeling  and  Wahl  (1997)  presented  the  results  of  SSI  analysis  of  the  pro¬ 
posed  North  Lock  Wall  at  McAlpine  Locks.  They  determined  that  the  downdrag 
force  was  significant,  and  that  it  could  be  substantially  affected  by  the  response 
of  the  interface  to  un  load-reload  cycles. 

Filz  and  Duncan  (1997)  presented  a  theory  to  quantify  the  downdrag  force  on 
the  back  of  nonmoving  retaining  walls.  Filz,  Duncan,  and  Ebeling  (1997)  pre¬ 
sented  a  simplified  method  for  incorporating  downdrag  forces  in  conventional 
analyses  of  nonmoving  retaining  walls.  They  also  observed  that  postconstruction 
settlement  causes  an  increase  in  the  downdrag  on  the  back  of  the  wall.  They  cited 
measurements  at  Eibach  Lock  in  Germany,  where  large  vertical  shear  forces  were 
persistent  for  10  years  under  repeated  filling  and  emptying  cycles  and  tempera¬ 
ture  fluctuations.  The  measured  Kv  remained  at  an  approximately  constant 
average  value  of  0.30.  These  vertical  shear  forces  cause  an  important  reduction  in 
the  lateral  earth  pressures  acting  on  the  wall. 

The  simplified  method  in  Filz,  Duncan,  and  Ebeling  (1997)  is  also  described 
in  Ebeling,  Pace,  and  Morrison  (1997).  This  method  is  shown  to  provide  an 
improvement  in  accuracy  over  conventional  analyses  of  nonmoving  walls  by 
producing  downdrag  values  that  are  closer  to  the  values  from  the  other  available 
sources  of  information:  (a)  theory,  (b)  finite  element  analyses,  (c)  data  from  pilot- 
scale  instrumented  retaining  wall  tests,  and  (d)  the  field  data  from  Eibach  Lock. 
While  more  accurate  than  conventional  analyses,  the  simplified  method  remains 
conservative  for  wall  stability  because  the  downdrag  values  from  the  simplified 


method  are  smaller  than  the  values  from  the  other  available  information. 
Furthermore,  the  evidence  from  the  instrumented  retaining  wall  tests  and  Eibach 
Lock  shows  that  downdrag  forces  either  remain  constant  or  increase  over  time 
after  backfilling.  The  simplified  method  is  restricted  to  backfill  materials  that  do 
not  exhibit  creep. 


1.3  Bonneville  Temporary  Tieback  Wall  Analysis 

1.3.1  Wall  description 

Strom  and  Ebeling  (2001)  used  the  Bonneville  temporary  tieback  wall  to 
illustrate  the  importance  of  the  NLFEM  in  the  analysis  and  design  of  projects. 
Portions  of  their  summary  are  contained  in  the  following  paragraphs.  The 
temporary  tieback  retaining  wall  is  approximately  134  m  (440  ft)  long.  The  wall 
is  constructed  by  slurry  trench  methods  in  6-m  (20-ft)-long  sections  similar  to 
those  shown  in  the  section  view  (horizontal)  of  Figure  1.4. 


Figure  1 .4.  Bonneville  navigation  lock,  temporary  tieback  wall-horizontal  section  (from  Strom  and 
Ebeling  2001) 


The  heights  of  the  panels  ranged  from  6  to  34  m  (20  to  1 1 0  ft).  Excavation 
sequencing  was  similar  to  that  shown  in  Figure  1.3.  Anchors  were  installed  in  a 
grid  pattern  of  approximately  3  m  (10  ft)  horizontal  by  3.4  m  (1 1  ft)  vertical. 

Each  tieback  anchor  was  composed  of  nineteen  15-mm  (0.6-in.)-diam  strands 
with  a  guaranteed  ultimate  strength  (GUTS)  of  1,862  MPa  (270  ksi).  Each  anchor 
was  prestressed  to  1 50  percent  of  its  design  load.  The  design  loads  are  approxi¬ 
mately  50  percent  of  the  anchor  ultimate  load  capacity.  Panel  6  was  the  focus  of 
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the  NLFEM  analysis.  A  section  view  (vertical)  of  Panel  6  is  shown  as  Figure  1 .5. 
The  tieback  anchor  loads  are  summarized  in  Table  1.1. 


Figure  1 .5.  Bonneville  navigation  lock,  temporary  tieback  wall-vertical  section  (from  Strom  and  Ebelinq 
2001)  a 


Table  1.1  =  ==  =— — =n 

Panel  6  Anchor  Loads  I 


Anchor 

Elevation 

Anchor 

Length 

Design  Load 

Prestress  Load 
at  150  percent 
Design  Load 

Lock-Off  Load 

m 

ft 

m 

ft 

kN 

1531 

KN 

159 

kN 

159 

25.6 

84 

27.1 

89 

1,249.9 

281 

1,874.9 

421.5 

1,209.9 

272 

22.3 

73 

24.1 

79 

1,249.9 

281 

1,874.9 

421.5 

1,298.8 

292 

18.9 

62 

20.7 

68 

1,249.9 

281 

1,874.9 

421.5 

1,289.9 

290 

15.5 

51 

15.8 

52 

1,592.5 

358 

2,388.7 

537.0 

1,583.6 

356 

Source: 

Strom  and  Ebeling  (2001). 

1.3.2  Overall  design  and  evaluation  process 

The  overall  design  and  evaluation  process  for  the  Bonneville  navigation  lock 
temporary  tieback  wall  involved  a  RIGID  analysis  and  a  WINKLER  analysis,  in 
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addition  to  the  NLFEM  described  in  this  chapter.  The  purposes  and  results  of  the 
RIGID  and  WINKLER  analyses  are  described  in  general  terms  and  illustrated 
using  Figure  1.6. 


Figure  1.6.  Overall  design  and  evaluation  process  (from  Strom  and  Ebeling  2001) 


1.3.3  RIGID  analysis 


A  simple  preliminary  analysis  using  the  RIGID  analysis  procedure  was 
®  employed  to  estimate  the  size  and  spacing  of  anchors.  The  RIGID  analysis  for 

Bonneville  was  based  on  a  composite  apparent  pressure  diagram,  as  shown  in 
Figure  1 .6a.  The  upper  rectangular-shaped  region  is  based  on  at-rest  pressures,  in 
accordance  with  Figure  29,  Chapter  3,  of  Design  Manual  7.2  (NAVFAC  1982). 
The  lower  triangular-shaped  region  is  also  based  on  at-rest  pressures  (Coulomb’s 
_  equation  with  a  factor  of  safety  of  1 .5  applied  to  the  shear  strength  of  the  soil). 

9  This  composite  diagram  assumed  that  wall  displacements  would  be  small. 
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thereby  keeping  earth  pressure  near  an  at-rest  state.  The  triangular  lower  portion 
of  the  diagram  was  considered  to  be  appropriate  since  the  upper  anchors  could 
possibly  lose  some  tension  during  the  life  of  the  wall,  thereby  leading  to  higher 
loads  at  the  base  of  the  wall.  Anchors  sizes  were  estimated  based  on  a  continuous 
beam  analysis  using  the  pressure  diagram  loading  of  Figure  1.6a.  There  was  very 
little  difference  between  the  anchor  loads  obtained  from  the  continuous  beam 
analysis  and  those  determined  using  the  tributary  area  method  (Munger,  Jones, 
and  Johnson  1990). 

The  wall  bending  moment  diagram  resulting  from  a  RIGID  analysis  is  illu¬ 
strated  in  Figure  1 ,6b.  Because  this  type  of  analysis  does  not  consider  the  wall 
displacement  that  occurs  during  prestressing  and  excavation,  it  will  underesti¬ 
mate  actual  bending  moment  demands  on  the  wall.  Therefore,  a  WINKLER 
analysis  was  performed  to  estimate  wall  bending  moments  and  to  confirm  anchor 
loads. 


1.3.4  WINKLER  analysis 

The  WINKLER  analysis  used  to  evaluate  the  Bonneville  navigation  lock 
temporary  tieback  wall  is  described  in  Munger,  Jones,  and  Johnson  (1990).  The 
moment  diagram  for  the  wall  was  determined  using  a  one-step  analysis  that 
modeled  the  final  excavation  condition  only.  The  moment  diagram  for  the 
WINKLER  analysis  is  illustrated  in  Figure  1.6c.  This  analysis,  because  it 
accounts  for  wall  displacements  during  anchor  prestressing  and  excavation,  does 
a  reasonably  good  job  of  estimating  moment  demands  on  the  wall.  A  RIGID 
analysis  that  considers  construction  staging  and  considers  support  displacements 
that  accumulate  during  each  stage  of  construction  (yielding  supports  analysis) 
could  also  provide  good  results. 

As  illustrated  in  Figure  1 ,6c,  the  shape  of  the  moment  diagram  obtained  by 
the  WINKLER  analysis  is  quite  different  from  that  obtained  by  the  standard 
RIGID  analysis.  Since  the  WINKLER  analysis  cannot  give  reliable  information 
with  respect  to  wall  displacements,  and  since  wall  displacements  and  displace¬ 
ments  of  the  ground  retained  by  the  wall  were  critical  to  performance,  it  was 
decided  to  perform  a  NLFEM  analysis.  The  NLFEM  analysis  investigated  wall 
displacements  and  bending  moment  demands  for  each  stage  of  construction. 


1.4  NLFEM  Analysis 

1.4.1  Ojectives 

The  purpose  of  the  NLFEM  analysis  was  threefold:  (1)  to  provide  a  means 
for  additional  confirmation  of  the  procedure  used  in  designing  the  wall,  (2)  to 
predict  potential  wall  performance  during  excavation  and  tieback  installation,  and 
(3)  to  assist  in  the  interpretation  of  instrumentation  results.  The  results  from  the 
study  depict  wall  behavior  in  terms  of  lateral  deflection,  bending  moments,  and 
earth  pressures  for  each  stage  of  construction. 


1.4.2  Description 

The  computer  program  SOILSTRUCT  (Ebeling,  Peters,  and  Clough  1992) 
was  used  for  the  analysis.  A  description  of  the  finite  element  analyses  can  be 
found  in  Mosher  and  Knowles  (1990).  SOILSTRUCT  is  designed  so  that  the 
actual  construction  process  can  be  simulated.  Simulation  of  the  actual  sequence 
of  construction  is  important  because  the  soil  stress  response  is  nonlinear  and 
stress  path  dependent.  SOILSTRUCT  provides  for  simulation  of  initial  stresses, 
fill  placement,  material  excavation,  dewatering,  and  placement  of  structural 
materials  in  a  series  of  incremental  load  steps.  Incremental  stresses  and  displace¬ 
ments  are  computed  after  each  load  step.  Table  1.2  lists  the  loading  steps  used  to 
model  the  sequence  for  the  wall  construction  and  lock  channel  excavation. 


Table  1.2 

Loading  Steps  in  SOILSTRUCT  Analysis 

Step 

Description 

1 

Construct  surcharge  to  pre-excavation  grade  (four  increments) 

2 

Excavate  for  railroad  relocation 

3 

Construct  slurry  trench  temporary  tieback  wall 

4 

Excavate  in  front  of  wall  to  el  78.5  ft  (23.9  m) 

5 

Install  upper  tieback  anchor  at  el  84  ft  (25.6  m)  and  prestress  to  1 50  percent  of  the 
design  load 

6 

Excavate  in  front  of  wall  to  el  67.5  ft  (20.6  m)  and  lock  off  upper  anchor  at  design  load 

7 

Install  second  tieback  anchor  at  el  73  ft  (22.3  m)  and  prestress  to  1 50  percentof  the 
design  load 

8 

Excavate  in  front  of  wall  to  el  56.5  ft  (17.2  m)  and  lock  off  second  anchor  at  design  load 

9 

Install  third  tieback  anchor  at  el  62  ft  (18.9  m)  and  prestress  to  150  percent  of  the  design 
load 

10 

Excavate  in  front  of  wall  to  el  45  ft  (13.7  m)  and  lock  off  third  anchor  at  design  load 

11 

Install  fourth  tieback  anchor  at  el  51  ft  (1 5.5  m)  and  prestress  to  1 50  percent  of  design 
load 

12 

Excavate  to  bottom  of  wall  at  el  39  ft  (11 .9  m)  and  lock  off  fourth  anchor  at  design  load 

I  Source:  Strom  and  Ebeling  (2001).  _ II 

1.4.3  Results 

The  results  for  each  stage  of  construction  were  studied,  from  the  in  situ  state 
and  wall  construction  through  the  excavation  and  tie  installation  procedure.  The 
results  for  each  stage  of  construction  are  illustrated  in  Figures  1.7a-1.7i  and 
described  below  in  terms  of  earth  pressures,  wall  bending  moments,  and 
displacements. 

1.4.3.1  Earth  pressures.  Lateral  earth  pressures  on  the  wall  for  each  stage 
of  construction  are  illustrated  in  Figures  1.7a-1.7i.  The  initial  pressure  on  the 
wall  is  approximately  50  percent  greater  than  at-rest  pressure.  This  increase  can 
be  attributed  to  overconsolidation  and  replacement  of  the  soil  by  a  concrete  wall. 
The  dotted  line  represents  at-rest  pressure  increased  by  50  percent.  The  earth 
pressure  distribution  changes  throughout  construction  as  a  result  of  excavation 
and  anchor  prestressing.  After  the  first  excavation  to  elevation  (el)  78.5  ft 
(23.9  m),  the  soil  behind  and  near  the  top  of  the  wall  is  in  an  active  state  as  a 
result  of  the  wall  moving  toward  the  excavation.  Farther  down  the  wall,  the 
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Figure  1.7.  Deflections,  moments,  and  earth  pressures  for  Bonneville  navigation  lock,  temporary  tieback 
wall-Panel  6  (from  Strom  and  Ebeling  2001 )  (Sheet  1  of  3) 


14 


Chapter  1  Introduction 


Chapter  1  Introduction 


Figure  1 .7.  (Sheet  3  of  3) 
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lateral  earth  pressure  is  greater  than  active,  but  less  than  the  initial  pressure  on 
the  wall.  The  analysis  shows  that,  with  prestressing  of  the  first  anchor,  earth 
pressures  increase  to  greater  than  initial  pressures  behind  the  upper  one  third  of 
the  wall.  Subsequent  excavations  and  anchor  prestressings  show  decreases  and 
increases,  respectively,  of  the  earth  pressures  along  the  wall.  Bulging  of  pressures 
around  the  anchors  appears  in  the  lower  one  half  of  the  wall.  Although  the  shape 
of  the  pressure  diagram  is  approximately  trapezoidal  during  the  initial  stages  of 
construction,  it  is  closer  to  triangular  for  the  later  stages  of  construction. 

1.4.3.2  Wall  bending  moments.  Wall  bending  moments  for  each  stage  of 
construction  are  illustrated  in  Figures  1.7a-1.7i.  Except  for  the  first  excavation, 
the  moment  diagram  for  the  wall  retains  the  same  general  form  throughout  con¬ 
struction.  The  region  immediately  behind  the  upper  anchor  experiences  negative 
bending  moment,  and  the  lower  region  below  the  upper  anchor  experiences  posi¬ 
tive  bending  moment.  The  maximum  moment  is  always  positive  and  varies 
during  construction  with  a  maximum  value  of  about  1 80  kN.m  (133  fit-kips).  The 
maximum  moment  develops  during  the  intermediate  stages  of  excavation  after 
the  second  anchor  is  prestressed. 

1.4.3.3  Wall  displacements.  Wall  displacements  for  each  stage  of  construc¬ 
tion  are  illustrated  in  Figures  1.7a-1.7i.  With  the  first  excavation  to  el  78.5  ft 
(23.9  m),  the  wall  moves  0.5  in.  (13  mm)  toward  the  excavation.  With  prestres¬ 
sing  of  the  upper  tieback  anchor,  the  wall  is  pulled  back  into  the  retained  soil, 
resulting  in  a  displacement  of  0.78  in.  (20  mm)  past  vertical  in  a  direction  away 
from  the  excavation.  In  subsequent  construction  steps,  there  is  little  change  in  the 
deflected  position  of  the  wall.  In  general,  the  wall  moves  into  the  soil  when  each 
anchor  is  prestressed  and  back  toward  vertical  with  each  excavation. 


1.4.4  Comparison  of  the  results  of  NLFEM,  WINKLER,  and 
RIGID  analyses 

The  analysis  of  the  Bonneville  navigation  lock  temporary  tieback  wall 
(presented  above)  demonstrates  the  significant  differences  in  the  results  of  the 
NLFEM,  WINKLER,  and  RIGID  approaches.  This  can  be  seen  for  instance  in 
the  assumed  apparent  earth  pressure  distribution  behind  the  wall,  and  the  moment 
diagrams  calculated  using  the  WINKLER  and  RIGID  procedures  (Figure  1 .6). 
These  diagrams  are  much  different  from  those  obtained  from  the  NLFEM  analy¬ 
sis  (Figure  1 .7),  which  shows  a  more  complicated  pattern  of  development  of 
earth  pressures  and  wall  bending  moments  with  construction  stage.  Prominent 
differences  are  in  nonuniform  development  of  earth  pressure  with  depth  and  the 
continuing  changes  in  the  magnitudes  of  the  earth  pressure  in  eveiy  construction 
stage.  In  contrast,  the  WINKLER  and  RIGID  procedures  assume  a  uniform  (until 
about  half  the  excavation  depth)  then  linearly  varying  earth  pressure  distribution 
with  depth.  There  is  no  soil  and  wall  interaction  effect  on  the  earth  pressure 
distribution.  The  moment  distributions  are  also  very  different  with  the  NLFEM 
analysis,  showing  reduction  in  bending  moment  after  installation  of  the  pre¬ 
stressed  anchors.  As  can  be  seen,  both  the  WINKLER  and  RIGID  procedures 
grossly  simplify  the  loads  acting  on  the  tieback  wall.  There  are  several  reasons 
for  the  differences  shown  in  Figures  1.6  and  1.7;  a  few  are  cited  below. 
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•  The  wall  is  assumed  to  be  stiff  in  the  RIGID  analysis  rather  than  flexible. 

•  Prestressing  to  achieve  stringent  displacement  control  objectives  results 
in  total  pressures  in  excess  of  at-rest  conditions. 

•  Construction  sequencing  has  a  major  impact  on  wall  performance. 

•  Pressure  distributions  are  approximately  triangular,  especially  during  the 
final  stages  of  construction. 

1 .5  Overview  of  SOILSTRUCT-ALPHA  Finite 
Element  Program 

Engineering  application  of  soil-structure  interaction  modeling  requires  a 
balance  between  modeling  realism  and  simplicity.  There  is  now  over  25  years  of 
experience  in  modeling  construction  procedures  by  the  finite  element  method, 
and  the  key  ingredients  in  engineering  applications  are  well  known.  One  of  the 
earliest  successful  applications  of  SSI  analysis  was  performed  by  Clough  and 
Duncan  (1969, 1971)  in  their  analysis  of  the  two  reinforced  concrete  U-frame 
locks  at  Port  Allen  and  Old  River.  These  two  locks  had  been  extensively  instru¬ 
mented.  Prior  to  Clough  and  Duncan’s  analyses,  the  instrumentation  data  had 
been  thought  to  be  unreliable  and  contrary  to  the  perceived  understanding  of  the 
behavior  of  locks  encountered  during  lock  operation.  Clough  and  Duncan’s  study 
showed  that  the  best  agreement  between  results  computed  using  the  finite  ele¬ 
ment  method  and  those  obtained  through  instrumentation  measurements  is 
obtained  when  the  actual  construction  process  is  simulated  as  closely  as  possible 
in  the  analysis. 

During  their  study,  Clough  and  Duncan  developed  what  is  referred  to  as  the 
backfill  placement  analysis  in  which  the  loads  exerted  by  the  backfill  on  the  lock 
wall  are  generated  automatically  during  simulated  placement  of  backfill  behind 
the  wall  (i.e.,  predetermined  earth  pressure  force  distributions  between  the  soil 
and  the  lock  are  not  specified).  This  requires  that  the  soil  backfill  and  foundation 
soil  strata  be  included  in  the  finite  element  mesh.  This  procedure  involves  the  use 
of  incremental  finite  element  analysis  with  nonlinear,  stress-dependent,  stress- 
strain  behavior  for  the  soil.  An  additional  requirement  is  that  interface  elements 
be  incorporated  within  the  finite  element  mesh  to  allow  for  relative  movement 
between  the  soil  and  structure. 

Clough  and  Duncan  developed  the  first  version  of  the  finite  element  program 
SOILSTRUCT,  which  implements  this  engineering  model ing/analysis  philoso¬ 
phy.  SOILSTRUCT  has  been  used  successfully  in  the  past  decade  on  several 
engineering  projects  supported  by  field  observations.  Since  1969,  several  ver¬ 
sions  of  SOILSTRUCT  have  been  developed  to  analyze  various  types  of  earth- 
retaining  structures  or  to  analyze  specific  problems  that  were  not  envisioned  at 
the  time  of  Clough  and  Duncan’s  original  development.  One  of  these  updated 
versions,  referred  to  as  SOILSTRUCT-ALPHA,  is  the  subject  of  this  report. 


SOILSTRUCT- ALPHA  (Ebeling,  Duncan,  and  Clough  1990)  is  a  special- 
purpose,  finite  element  program  for  two-dimensional  (2-D),  plane  strain  analysis 
of  SSI  problems.  SOILSTRUCT  calculates  displacements  and  stresses  resulting 
from  incremental  construction,  backfilling,  excavation,  dewatering,  rising  water 
table,  and/or  load  application.  Nonlinear,  stress-path  dependent,  stress-strain 
behavior  of  the  backfill  was  approximated  in  the  finite  element  analysis  using  the 
tangent  modulus  method.  In  the  tangent  modulus  method,  new  values  of  tangent 
moduli  are  assigned  to  each  soil  element  at  each  increment  of  loading  (i.e., 
dewatering,  lock  construction,  and  backfilling)  or  unloading  (i.e.,  excavation, 
rising  water  table).  The  modulus  values  assigned  to  each  element  are  adjusted  in 
accordance  with  their  stresses  to  simulate  nonlinear  behavior. 

Three  types  of  finite  elements  are  used  to  represent  the  behavior  of  different 
materials: 

1.  Two-dimensional  continua  elements 

A  2-D,  subparametric,  quadrilateral  element  (QM5)  is  used  to  represent  the 
soil  and  most  structural  materials.  Structural  supports,  such  as  the  struts  or 
tieback  components  of  an  excavation  support  system,  are  typically  modeled 
as  a  spring  support  using  bar  elements.  However,  2-D  elements  have  been 
used  to  model  these  supports  in  some  cases. 

2.  Interface  elements 

SOILSTRUCT- ALPHA  has  the  ability  to  model  the  interface  region  between 
the  soil  backfill  and  the  structure  using  interface  elements.  This  important 
feature  allows  for  the  movement  of  the  softer  continua  elements  used  in 
modeling  the  backfill  relative  to  (the  movement  of)  the  stiffer  continua 
elements  used  in  modeling  the  structure.  This  element  is  defined  by  four 
nodes,  with  each  of  the  two  pairs  of  nodes  having  the  same  coordinates;  thus 
this  element  has  no  thickness. 

3.  One-dimensional  bar  elements 

To  model  the  behavior  of  a  variety  of  structural  systems,  1-D,  two-node,  bar 
or  spring  elements  are  used.  This  includes  the  modeling  of  structural  supports 
such  as  braces  or  tiebacks  or  the  modeling  of  reinforcement  placed  within  a 
soil  backfill. 

SOILSTRUCT  was  expanded  during  the  U.S.  Corps  of  Engineers’  Repair, 
Evaluation,  Maintenance,  and  Rehabilitation  Research  Program  to  model  the  loss 
of  contact  between  the  base  of  a  wall  (a  lock,  in  this  case)  and  its  rock  foundation 
using  a  procedure  called  the  ALPHA  method  (Ebeling,  Duncan,  and  Clough 
1990;  Ebeling  et  al.  1992).  Further  details  on  the  ALPHA  method  are  given  in 
Ebeling  et  al.  (1992,  pages  64-70).  The  ALPHA  method  was  extended  to  soil 
elements  by  Regaldo,  Duncan,  and  Clough  (1992)  to  reduce  numerical 
inaccuracies  in  soil  elements  that  are  at  or  near  failure. 

The  continua  elements  used  to  model  the  soil  and  the  soil-to-structure  inter¬ 
face  elements  that  may  have  failed  in  shear  at  one  stage  of  loading  have  the 
ability  to  recover  their  shear  stiffness  and  shearing  resistance  as  a  result  of  an 
increase  in  confining  pressures  at  some  later  stage  of  loading  in  this  version  of 
SOILSTRUCT-ALPHA.  Several  other  improvements  have  been  made  to  the 
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material  models,  including  the  new  Gomez,  Filz,  and  Ebeling  (2000a,  2000b) 
extended  load/unload/reload  hyperbolic  model  for  interfaces,  and  to  the  numeri¬ 
cal  procedures  implemented  within  SOILSTRUCT-ALPHA  based  on  experience 
gained  at  the  U.S.  Army  Engineer  Research  and  Development  Center  in  conduct¬ 
ing  SSI  analyses  of  different  types  of  structures. 


1 .6  Need  for  Further  Improvements  in 
SOILSTRUCT-ALPHA 

1.6.1  Improvements  in  excavation  procedure 

The  Corps  uses  SOILSTRUCT-ALPHA  to  perform  SSI  analyses  of  multi- 
anchored  or  tieback  retaining  walls.  Permanent  multi-anchored  walls  have  been 
used  as  guide  walls  and  approach  walls  on  navigation  projects,  and  as  retaining 
walls  on  highway  and  railroad  protection  and  relocation  projects  (Mosher  and 
Knowles  1990).  Multi-anchored  walls  are  constructed  by  first  installing  soldier 
elements  for  entire  wall  panels,  and  then  beginning  excavation.  Anchors  or 
tiebacks  are  installed  as  excavation  proceeds.  A  critical  aspect  of  SSI  analysis  of 
such  structures  is  modeling  the  excavation  process. 

The  excavation  algorithm  in  SOILSTRUCT  was  developed  in  the  late  1960s 
by  Clough  and  Duncan  (1969),  and  it  has  not  been  updated  since  that  time,  even 
though  improved  methods  have  been  developed  by  Ghaboussi  and  Pecknold 
(1984),  Borja,  Lee,  and  Seed  (1989),  and  others.  The  Clough  and  Duncan  algo¬ 
rithm  is  based  on  extrapolation  of  stresses  from  nearby  elements  to  the  excava¬ 
tion  boundary  and  integration  of  the  extrapolated  stresses  to  produce  excavation 
unloading  forces  at  element  nodes.  One  of  the  weaknesses  of  the  Clough  and 
Duncan  excavation  algorithm  is  that  the  forces  applied  along  the  excavation 
boundary  are  not  always  consistent  with  the  stresses  in  the  excavated  elements. 
Force  equilibrium  is  not  necessarily  satisfied,  and  the  resulting  deformations  can 
be  either  larger  or  smaller  than  they  should  be.  The  lack  of  force  equilibrium  can 
result  in  a  distribution  of  excavation  unloading  force  in  which  either  too  much  or 
too  little  force  is  applied  to  the  underlying  soil  and,  correspondingly,  either  too 
little  or  too  much  force  is  applied  to  the  structure.  Moreover,  it  has  been  shown 
that  the  Clough  and  Duncan  procedure  may  yield  results  that  are  dependent  on 
the  number  of  excavation  stages  used  in  simulations  for  linear  elastic  materials. 
For  such  materials,  it  has  been  shown  that  the  results  of  the  excavation  should  be 
independent  of  the  number  of  steps  and  the  sequence  used  to  simulate  the 
excavation. 


1.6.2  Accounting  for  soil-structure  interface  behavior  in  the 
excavation  procedure 

A  substantial  amount  of  research  has  been  performed  in  recent  years  on  the 
interaction  of  soils  and  structures  along  their  interface  (Ebeling  et  al.  1993; 
Ebeling  and  Mosher  1996;  Ebeling  and  Wahl  1997;  Ebeling,  Pace,  and  Morrison 
1997;  Ebeling,  Peters,  and  Mosher  1997).  These  studies  showed  that  the  behavior 
of  the  soil-structure  interface  has  a  significant  influence  on  the  magnitudes  of  the 
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loads  acting  against  earth-retaining  structures  such  as  lock  walls  and  soil  rein¬ 
forcements  such  as  slurry  walls.  They  also  illustrated  that  the  pre-  and  post¬ 
construction  stress  paths  followed  by  interface  elements  are  complex,  often 
involving  simultaneous  changes  in  normal  and  shear  stresses,  as  well  as 
unloading-reloading  due  to  postconstruction  rise  of  the  groundwater  level. 

Gomez,  Filz,  and  Ebeling  (2000a)  developed  an  extended  hyperbolic  model 
for  interfaces  and  implemented  it  into  the  finite  element  program  SOILSTRUCT- 
ALPHA.  The  model  is  based  on  the  Clough  and  Duncan  (1971)  hyperbolic 
formulation,  which  was  extended  to  model  a  variety  of  stress  paths.  Gomez,  Filz, 
and  Ebeling  (2000a)  performed  a  series  of  interface  tests  between  uniform  fine 
sands  and  concrete.  Some  of  these  tests  followed  complex  stress  paths  that 
included  unloading-reloading  and  simultaneous  changes  in  normal  and  shear 
stresses. 

Gomez,  Filz,  and  Ebeling  (2000a)  also  carried  out  a  pilot-scale  lock  wall 
simulation  that  modeled  placement  and  compaction  of  the  backfill,  surcharge 
application,  and  changes  in  the  elevation  of  the  water  table  behind  the  wall.  By 
comparing  model  predictions  to  interface  test  results,  and  the  results  of 
SOILSTRUCT-ALPHA  analyses  to  measurements  from  the  lock  wall  simulation, 
these  investigators  concluded  that  the  extended  load/unload/reload  hyperbolic 
model  could  provide  accurate  estimates  of  the  response  of  backfill-to-lock  wall 
interfaces. 

Important  similarities  exist  between  the  types  of  loading  that  occur  at  struc- 
ture-to-soil  interfaces  in  both  multi-anchored  systems  and  lock  walls.  Therefore, 
it  is  possible  that  the  Gomez-Filz-Ebeling  model  for  interfaces  could  also  be  used 
for  SSI  analyses  of  multi-anchored  systems.  However,  because  the  model  was 
developed  based  on  the  results  of  interface  tests  performed  using  uniform  fine 
sands,  additional  testing  was  required  to  validate  model  performance  for  coarser 
soils. 

Gomez,  Filz,  and  Ebeling  (2000b)  performed  a  series  of  virgin  shear  tests 
under  constant  stress  at  the  interface  between  a  coarse  sand  and  concrete.  The 
results  of  these  tests  were  used  to  determine  the  hyperbolic  parameter  values  of 
the  interface  following  the  recommendations  given  by  Gomez,  Filz,  and  Ebeling 
(2000a).  An  interface  test  was  performed  following  a  complex  stress  path  that 
included  unloading-reloading  as  well  as  simultaneous  changes  in  shear  and  nor¬ 
mal  stresses.  The  interface  response  measured  during  this  test  was  compared  with 
the  response  calculated  using  the  extended  hyperbolic  model.  It  was  found  that 
the  Gomez-Filz-Ebeling  interface  model  provided  accurate  estimates  of  the 
response  of  this  type  of  interface.  Therefore,  it  can  be  concluded  that  the 
extended  hyperbolic  model  can  be  used  for  prediction  of  the  response  of  inter¬ 
faces  between  concrete  and  a  variety  of  granular  soils.  The  hyperbolic  parameter 
values  of  the  interface  tested  also  add  to  the  database  of  interface  properties 
available  in  the  literature.  The  extended  hyperbolic  model,  together  with  the 
interface  data  that  have  been  generated,  provides  a  useful  tool  for  analyzing 
multi-anchored  retaining  systems  and  other  Corps  of  Engineers  structures. 

The  excavation  algorithm  in  SOILSTRUCT-ALPHA  was  not  developed  with 
consideration  of  interface  elements  between  the  excavated  soil  and  the  retaining 
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structure.  Such  interface  elements  are  necessary  for  realistic  modeling  of  earth- 
retaining  structures  because  the  soil  has  a  tendency  to  slip  at  its  contact  with  the 
retaining  structure,  and  interface  elements  allow  slip  to  occur  in  finite  element 
analyses.  Consequently,  it  is  important  to  upgrade  SOILSTRUCT-ALPHA  with 
an  excavation  procedure  that  correctly  accounts  for  the  presences  of  interface 
elements. 


1 .7  Study  Objectives 

The  purpose  of  this  report  is  to  provide  specific  recommendations  for 
improving  the  excavation  algorithm  in  SOILSTRUCT,  and  to  provide  the  infor¬ 
mation  necessaiy  for  implementing  these  recommendations.  The  results  will 
permit  reliable  and  efficient  improvement  of  the  excavation  algorithm  in 
SOILSTRUCT-ALPHA.  The  objectives  of  the  study  described  in  this  report 
are  to 

a.  Complete  a  comprehensive  literature  review  of  numerical  modeling  of 
excavation. 

b.  Determine  the  suitability  of  existing  excavation  algorithms  for  their  use 
in  SOILSTRUCT. 

c.  Provide  the  information  necessaiy  to  implement  an  improved  excavation 
algorithm  in  SOILSTRUCT. 

1.8  Report  Organization 

The  report  is  organized  in  five  chapters  and  one  appendix.  Chapter  1 
provides  a  general  background  on  the  geotechnical  issues  addressed  in  the  study 
and  the  motivation  for  the  work  carried  out  in  the  study. 

Chapter  2  presents  a  detailed  discussion  of  the  general  numerical  procedure 
to  simulate  excavation.  The  chapter  focuses  on  the  Clough  and  Duncan  pro¬ 
cedure,  and  the  use  of  the  procedure  is  illustrated  via  an  example  problem. 

Chapter  3  gives  a  review  of  the  other  available  numerical  procedures  for 
simulating  excavation.  Four  major  types  of  numerical  procedures  are  identified, 
and  their  main  features  are  discussed.  A  detailed  presentation  of  the  method  of 
force  residuals  is  given  by  analyzing  the  example  problem  used  in  Chapter  2. 

A  comprehensive  assessment  and  comparison  of  the  available  numerical 
models  for  excavation  is  given  in  Chapter  4.  The  models  are  compared  in  terms 
of  practical  and  theoretical  advantages  and  disadvantages.  The  Clough  and 
Duncan  procedure  and  the  procedure  based  on  force  residuals  are  assessed  with 
regard  to  how  and  whether  they  satisfy  the  requirements  of  force  equilibrium, 
uniqueness  of  solution,  and  accuracy  and  convergence  of  solution.  Chapter  4  also 
gives  a  detailed  discussion  of  the  issues  related  to  the  implementation  of  the 
method  of  force  residuals  in  SOILSTRUCT-ALPHA. 
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Chapter  5  summarizes  the  conclusions  from  the  study  and  gives  recommen¬ 
dations  for  updating  and  modifying  the  excavation  algorithm  used  in 
SOILSTRUCT-ALPHA. 

Finally,  Appendix  A  gives  a  brief  discussion  on  the  determination  of  external 
surface  tractions  and  internal  element  forces  using  isoparametric  finite  elements. 
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2  Numerical  Simulation  of 
Excavation 


Typical  examples  of  excavation  in  soils  were  shown  in  Chapter  1 .  In  finite 
element  analysis,  these  excavations  are  simulated  following  the  basic  steps 
shown  in  Figure  2.1 .  The  particular  example  used  in  the  illustration  is  the  excava¬ 
tion  stage  during  construction  of  a  slurry  wall  along  a  river  embankment  (see 
Figures  1 ,3b-c).  The  portion  of  the  soil  to  be  excavated  is  shown  as  a  shaded 
area.  Prior  to  excavation,  the  soil  to  be  excavated  exerts  stresses  along  the 
boundary  between  the  soil  to  be  excavated  and  the  remaining  soil  (Figure  2.1a). 
As  far  as  the  remaining  soil  is  concerned,  the  situation  is  the  same  as  in  Fig¬ 
ure  2.1b  where  the  excavated  soil  is  simply  replaced  by  the  stresses  (or  tractions) 
across  the  face  of  the  excavated  surface.  The  soil  that  will  remain  after  the  exca¬ 
vation  would  undergo  no  displacement  or  change  in  stress  if  the  soil  to  be  exca¬ 
vated  were  replaced  by  the  boundary  stresses.  Since  the  systems  shown  in 
Figures  2.1a-b  are  equivalent,  excavation  simply  involves  removal  of  the  bound¬ 
ary  stresses  from  the  remaining  soil  (Figure  2.1c).  A  stress-free  surface  would 
then  result,  and  the  displacements  and  stresses  due  to  excavation  would  be  pro¬ 
duced  by  the  removal  of  this  load.  In  summary,  simulation  of  excavation  involves 
the  following  steps: 

a.  Find  the  tractions  or  boundary  stresses  transmitted  to  the  remaining  soil 
by  the  soil  that  will  be  excavated. 

b.  Remove  the  stiffness  of  the  excavated  region  from  the  stiffness  of  the 
whole  region. 

c.  Apply  to  the  remaining  soil  the  tractions  or  boundary  stresses  with 
magnitudes  equal  to  those  determined  in  Step  1  and  opposite  in  sign. 

d.  Add  the  incremental  displacements,  strains,  and  stresses  of  Step  3  to  the 
condition  before  excavation. 

In  finite  element  simulation,  stresses  induced  by  external  loads  and  self-weight 
are  converted  to  nodal  forces. 
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Figure  2.1.  Steps  in  simulation  of  excavation 
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2.1  Clough  and  Duncan  Procedure 

One  of  first  finite  elements  procedures  for  simulating  excavation  in  soils 
was  developed  by  Clough  and  Duncan  (1969).  This  is  the  simulation  procedure 
used  in  the  finite  element  code  SOILSTRUCT  and  the  latest  version  of 
SOILSTRUCT-ALPHA.  Subsequent  applications  showed  that  the  Clough  and 
Duncan  (1969)  procedure  suffers  from  certain  limitations,  which  could  have 
significant  influence  on  the  simulated  results.  Following  the  work  of  Clough  and 
Duncan  (1969),  several  other  numerical  procedures  for  modeling  excavation  have 
been  developed  that  improved  upon  and  removed  the  deficiencies  in  the  Clough 
and  Duncan  procedure.  This  section  gives  a  detailed  review  of  the  Clough  and 
Duncan  procedure.  The  review  is  presented  to  provide  a  starting  point  for  asses¬ 
sing  alternative  methods  for  modeling  excavation  and  determining  how  such 
alternative  methods  can  be  implemented  in  a  modified  and  updated  version  of 
SOILSTRUCT-ALPHA. 

Similar  to  many  excavation  procedures,  the  technique  of  Clough  and  Duncan 
is  based  on  Figure  2.1  and  consists  of  determining  the  stresses  at  the  excavation 
boundary,  calculating  the  equivalent  nodal  forces  from  the  boundary  stresses,  and 
applying  equal  but  opposite  nodal  loads  to  the  finite  element  mesh.  At  the  same 
time,  the  contributions  of  the  excavated  elements  are  taken  out  of  the  global 
stiffness  matrix.  This  can  be  done  by  assigning  minimal  stiffness  values  to  the 
excavated  elements.  The  application  of  equal  but  opposite  loads  at  the  nodes 
along  the  excavation  boundary  ensures  that  the  excavation  boundary  becomes 
stress  free  following  the  removal  of  the  excavated  elements. 

To  calculate  the  nodal  loads  on  the  excavation  boundaries,  the  Clough  and 
Duncan  procedure  uses  the  normal  and  shear  stresses  along  the  surface  to  be 
exposed  by  the  excavation.  Consider  the  case  of  a  four-noded  finite  element  with 
a  linear  distribution  of  normal  and  shear  stresses  along  the  excavation  edges  of 
the  element,  as  illustrated  in  Figure  2.2.  The  normal  stresses  will  produce  nodal 
forces  perpendicular  to  the  excavation  edge,  while  the  shear  stresses  will  produce 
nodal  forces  parallel  to  the  excavation  edge.  Addition  of  all  forces  in  the  same  x 
and  y  directions  gives  the  total  excavation  forces  Fx  and  Fy  at  each  node. 

The  nodal  forces  for  the  case  of  uniform  stress  distribution  along  the  excava¬ 
tion  edge  will  be  equal  to  one  half  the  stress  multiplied  by  the  length  of  the  exca¬ 
vated  edge  for  both  nodes  (Figure  2.3).  For  quadrilateral  elements  with  linear 
variation  of  stress  along  the  edges,  the  nodal  forces  are  calculated  as  shown  in 
Figure  2.3.  The  nodal  loads  are  obtained  from  simple  static  equilibrium  and  by 
integrating  the  area  of  the  load  over  the  length  of  side  of  the  element.  Appen¬ 
dix  A  presents  a  derivation  of  the  equivalent  nodal  loads  given  in  Figure  2.3. 

In  the  finite  element  method,  the  stresses  are  generally  determined  only 
within  the  element  (at  the  center  of  the  elements  or  at  the  Gaussian  quadrature 
points)  and  not  along  element  edges.  Since  excavation  boundaries  pass  between 
elements,  several  procedures  were  developed  to  determine  the  stresses  on  the 
excavation  boundaries  from  the  internal  stresses  of  the  elements  adjacent  to  the 
excavation  boundary.  Dunlop,  Duncan,  and  Seed  (1968)  determined  the  stresses 
on  the  boundaries  by  averaging  the  stresses  in  pairs  of  elements  adjacent  to  the 
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Figure  2.2.  Normal  and  shear  stresses  and  equivalent  nodal  forces  along  the  excavated  sides  of  an 
element 
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b.  Linear  distribution 


Figure  2.3.  Equivalent  nodal  forces  from  surface  stresses 


excavation  boundary.  Equivalent  nodal  forces  were  then  calculated  assuming  that 
the  average  stresses  are  constant  along  the  excavation  boundary.  In  the  case 
shown  in  Figure  2.4,  the  stresses  along  the  excavation  boundary  K-L  are  obtained 
from  the  average  of  the  stresses  in  Elements  1  and  2,  and  for  the  boundary  J-K 
the  stresses  are  calculated  from  Elements  1  and  4.  This  method  was  shown  to  be 
accurate  only  when  the  elements  are  of  equal  size.  Chang  (1969)  improved  on 
this  method  by  using  only  the  stresses  above  the  boundary  and  correcting  the 
nodal  load  according  to  an  assumed  gravity  stress  gradient  within  the  elements. 

Clough  and  Duncan  generalized  the  techniques  of  Dunlop,  Duncan,  and  Seed 
(1968)  and  Chang  (1969)  by  interpolating  the  stresses  along  the  excavation 
boundary  from  the  stresses  in  the  element  to  be  excavated  and  the  stresses  in  the 
three  adjacent  elements.  This  is  the  procedure  that  is  currently  used  in 
SOILSTRUCT-ALPHA.  Figure  2.4  illustrates  the  Clough  and  Duncan  procedure. 
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Figure  2.4.  Interpolation  procedures  used  to  determine  the  stresses  along  the 
excavation  boundaries.  For  the  Dunlop,  Duncan,  and  Seed  (1968) 
procedure,  stresses  along  K-L  are  calculated  from  the  stresses  of 
Elements  1  and  2,  and  stresses  along  J-K  are  calculated  from  the 
stresses  of  Elements  1  and  4.  For  the  Clough  and  Duncan  (1969) 
procedure,  stresses  at  the  nodal  points  J,  K,  and  L  are  interpolated 
from  the  stresses  in  the  centers  of  the  Elements  1,2,3,  and  4. 
Elements  2,  3,  and  4  adjacent  to  the  excavated  Element  1  are  called 
interpolation  elements 


To  interpolate  the  stresses  along  the  excavation  boundary,  it  is  assumed  that  the 
stresses  vary  bilinearly  in  the  region  near  the  excavation  boundary.  An  inter¬ 
polation  function  of  the  following  form  is  used  to  express  the  variation  in  stress: 

o  =  al+a2x  +  a3y  +  a4xy  (2.1) 


where  x  and  y  are  the  coordinates  of  a  nodal  point  and  o  is  any  of  the  stress 
components  ox,  cy,  and  <5xy .  The  interpolation  coefficients  a, ,  a2 ,  a3,  and  a4 

are  determined  from  four  known  values  of  the  stress  cr  at  four  locations  given  by 
the  coordinates  x  andy.  As  mentioned  above,  these  four  locations  are  the  centers 
of  the  excavated  element  and  the  three  adjacent  elements  (the  three  adjacent 
elements  are  called  the  interpolation  elements).  Numbering  the  elements  1  to  4, 
the  stresses  o(l),  a(2),  c(3),  and  o(4)  in  each  element  can  be  expressed  in  terms 
of  the  coordinates  of  the  center  of  the  elements  (jci,  yO  . . .  (x4,  y4)  and  the 
interpolation  coefficients  ax,  a, ,  a3 ,  and  a4  as 

o(l)  =  ax  +  a2x]  +  a3yx  +  a4xxyx 

0(2)  =  a,  +  a2x2  +  a3y2  +  a4x2y2  ^  2 

o(3)  =  a,  +  a2x3  +  a3y3  +  a4x3y3 
o(4)  =  a,  +  a2x4  +  a3y4  +  a4x4y4 
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These  equations  can  be  written  in  matrix  form  as 

{<*},=  MM 


(2.3) 


where  {<?},=  {a(l)  a(2)  a(3)  o(4)}r ,  {a}  =  {a,  a2  a3  a4}r,and 


M  = 


1  *  y>  w 

1  *2  y2  x2y2 

1  y3  w 

.1  xa  y4  x4y4 


(2.4) 


The  unknown  interpolation  coefficients  {a}  =  {a,  a2  a,  a4)'  can  be 
obtained  from  Equation  2.2  as 


{a}  =  [/n]  '  {<j} 


(2.5) 


The  values  of  the  interpolation  coefficients  can  now  be  used  to  extrapolate 
the  stresses  at  each  of  the  nodes  of  the  element  to  be  excavated.  Denoting  the 
four  nodes  of  the  element  to  be  excavated  as  I,  J,  K,  and  L,  the  stresses  at  the 

nodes  {c}n  =  {a;  cy  aK  oL  }r  can  be  determined  as 

H=MM  (2.6) 


where  [n]  is  the  coordinate  matrix  for  the  nodes  I,  J,  K,  and  L: 


X, 

y, 

x,y, 

XJ 

yj 

xjyj 

XK 

yK 

Wk 

XL 

y,. 

xi.y,. 

(2.7) 


Finally,  the  nodal  forces  at  the  excavation  boundary  nodes  can  be  calculated 
from  the  nodal  stresses  assuming  the  stresses  vary  linearly  along  the  sides  of  the 
element.  The  calculations  of  the  nodal  forces  from  the  surface  tractions  are  the 
same  as  those  shown  in  Figure  2.3. 

A  three-dimensional  presentation  of  the  Clough  and  Duncan  interpolation 
procedure  is  shown  in  Figure  2.5.  Essentially,  the  procedure  fits  a  bilinear  distri¬ 
bution  of  stresses  from  the  interpolation  elements  and  uses  the  interpolation 
parameters  in  finding  the  stresses  at  the  nodes  of  the  excavated  element.  This 
means  that  the  interpolated  surface  is  extended  from  the  interpolation  elements  to 
the  region  occupied  by  the  excavated  elements. 
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The  same  surface  is  applied 
to  determine  the  stresses  at 
the  nodal  points  of  the 
excavated  elements. 
Boundary  stresses  and 
excavation  forces  are  then 
calculated  assuming  linear 
distribution  of  stress  between 
the  nodes. 


•  center  of  interpolation 
element 

♦  nodes  of  excavated 
element 


▲ 

Surface  interpolated  from 

stresses  in  the  centers  of  the 

four  interpolation  elements: 

_ — ^ 

o  =  ax  +  a2x  +  a3y  +  aAxy 

Figure  2.5.  Three-dimensional  representation  of  the  Clough  and  Duncan 
interpolation  procedure 


The  Clough  and  Duncan  procedure  is  illustrated  in  the  flowchart  shown  in 
Figure  2.6. 


2.2  Example  Calculation 

To  illustrate  the  Clough  and  Duncan  procedure,  the  simple  problem  shown  in 
Figure  2.7,  which  consists  of  an  excavation  in  level  ground,  will  be  used.  For 
convenience,  constant  stress  rectangular  elements  will  be  used.  The  same  prob¬ 
lem  will  also  be  used  in  Section  3.3  to  illustrate  methods  based  on  the  calculation 
of  excavation  forces  using  internal  element  stresses. 


2.2.1  Existing  ground  condition  (pre-excavation) 

The  ground  consists  of  a  2-m  horizontal  layer  of  homogenous,  normally 
consolidated  loose  sand  underlain  by  stiff  rock.  For  convenience,  the  unit  weight 
of  the  soil  is  assumed  to  be  y  =  20  kN/m3.  The  initial  in  situ  vertical  and  hori¬ 
zontal  stresses  are  due  to  the  self- weight  of  the  soil.  The  in  situ  vertical  stress  oy 

varies  linearly  with  depth  y: 

Gy=yy  (2-8) 
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Figure  2.6.  Flowchart  for  excavation  modeling  using  the  Clough  and  Duncan 
procedure 

Assuming  an  infinite  horizontal  extent,  uniaxial  oedometric  conditions  apply. 
Consequently,  the  horizontal  <yx  is  related  to  the  vertical  stress  according  to 

(2.9) 


where  Ka  is  the  at-rest  ratio  of  horizontal  to  vertical  stress.  This  ratio  is  estimated 
from  Jaky’s  equation  for  normally  consolidated  deposits  as 
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Figure  2.7.  Initial  in  situ  stresses  beneath  a  level  ground  surface 


K0  =  1  —  sintj) 


(2.10) 


Assuming  the  friction  angle  <|>  is  equal  to  30  deg  gives  Ka  =0.5.  Using  y  = 
20  kN/m3  and  Ka  =  0.5  gives  the  in  situ  vertical  and  horizontal  stress 
distributions  shown  in  Figure  2.7. 


2.2.2  Finite  element  analysis 

^  Step  1:  Gravity  turn-on  analysis.  The  first  step  is  the  simulation  of  the 

ground  conditions  of  the  site  before  the  start  of  the  construction  operations.  In 
finite  element  programs  such  as  SOILSRUCT  and  SOILSTRUCT-ALPHA,  the 
standard  procedure  involves  representing  the  initial  conditions  by  an  initial  state 
of  stress  in  each  element  in  which  the  vertical  stress  is  equal  to  the  overburden 
•  stress,  and  the  horizontal  stress  is  a  fraction  of  the  vertical  stress  given  by  Equa¬ 

tion  2.9.  Assuming  level  ground,  the  initial  shear  stresses  on  the  horizontal  and 
vertical  planes  are  assumed  to  be  zero.  Also,  the  displacements  and  strains  from 
the  initial  stresses  are  set  to  zero  once  the  gravity  turn-on  analysis  is  completed. 
However,  the  stresses  calculated  in  subsequent  incremental  loadings  are  cumula¬ 
tively  superimposed  on  the  initial  stresses.  It  is  important  to  account  for  the  initial 
^  stresses,  particularly  in  cases  where  stress-dependent  deformation  moduli  are 

used.  This  gravity  turn-on  procedure  does  not  explicitly  model  the  influence  of 
stress  history  on  material  behavior. 

A  simplified  finite  element  model  of  the  level  ground  shown  in  Figure  2.7  is 
_  given  in  Figure  2.8.  The  finite  element  model  consists  of  two  layers  of  constant 

9  stress  elements  and  a  total  of  six  elements.  The  calculation  involves  the  determin¬ 

ation  of  the  nodal  loads  at  the  excavation  boundaries  due  to  the  removal  of  Ele¬ 
ment  2.  To  simulate  an  infinite  horizontal  extent,  no  lateral  displacements  are 
allowed  along  the  two  sides.  Also  since  the  layer  is  underlain  by  rock,  the  vertical 
displacement  of  the  base  of  the  model  is  assumed  to  be  zero.  The  size  of  all 
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Figure  2.8.  Excavation  of  a  single  element 


elements  is  1  by  1  m2,  and  a  unit  thickness  is  assumed.  Elastic  properties  are 
assumed  (see  below). 

Because  loading  in  finite  element  analysis  can  only  be  represented  by  the 
application  of  nodal  forces,  the  self-weights  of  the  elements  have  to  be  converted 
to  vertical  nodal  forces.  Again,  assuming  a  unit  thickness  of  the  model,  the  self¬ 
weight  of  each  element  W  shown  in  Figure  2.8  is  equal  to 

W  =  yL2=( 20  kN/m3)  (1  m)(lmXlm)  =  20  kN  (2.11) 


For  a  rectangular  mesh,  the  equivalent  nodal  loads  from  the  element  self-weights 
are  obtained  by  adding  one  quarter  of  the  weight  of  each  element  and  summing 
the  weights  of  all  elements  connected  to  a  node  (Cook,  Malkus,  and  Plesha 
2002).  For  the  problem  shown  in  Figures  2.7  and  2.8,  the  self-weight  induces  the 
vertical  nodal  load  vectors  given  below  (Figure  2.9).  For  example,  the  vertical 
load  of  10  kN  for  node  2  comes  from  one  quarter  the  weight  of  Element  1  (equal 
to  5  kN)  and  one  quarter  the  weight  of  Element  2  (equal  to  5  kN).  The  summa¬ 
tion  of  element  weights  is  carried  out  for  all  the  nodes. 

Step  2:  Results  of  gravity  turn-on  analysis.  Having  prescribed  the  vertical 
loads  due  to  the  element  self-weights,  a  finite  element  analysis  can  now  be 
carried  out  to  calculate  the  initial  stresses  in  the  elements.  Normally,  this  will 
also  require  that  the  material  parameters  be  prescribed.  For  an  elastic  material,  a 
Young’s  modulus  E  and  a  Poisson’s  ratio  v  are  required.  For  the  situation  shown 
in  Figure  2.9,  the  analysis  of  the  element  vertical  stresses  is  trivial  and  can  be 
done  solely  based  on  vertical  equilibrium.  The  Young’s  modulus  is  assumed  to 
be  1 0  MPa. 
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Figure  2.9.  Nodal  loads  from  self-weight  (forces  shown  are  in  kilonewtons) 

The  horizontal  stresses  can  be  obtained  from  an  elastic  stress-strain  relation 
assuming  the  horizontal  strain  is  equal  to  zero.  The  horizontal  stresses  from 
this  elastic  stress-strain  solution  follow  Equation  2.9,  where  K0  is  related  to  the 
Poisson’s  ratio  v  of  the  material  by  the  following  relation: 

Ka= —  (2-12) 

°  1 — v 


Thus,  to  achieve  a  Ka  value  consistent  with  the  value  obtained  from 
Equation  2.10,  the  Poisson’s  ratio  v  can  be  solved  from  Equation  2.9  as 


v  = 


K„ 

l  +  K0 


(2.13) 


This  gives  v  =  0.33  for  Ka  =  0.5.  Based  on  vertical  equilibrium  and  the  solution 

of  the  horizontal  elastic  stress-strain  relation,  the  vertical  and  horizontal  stresses 
at  the  centers  of  the  elements  (due  to  the  self- weight  vertical  loads  shown  in 
Figure  2.9)  are  given  in  Table  2.1  below. 


I  Table  2.1 

Calculated  Element  Stresses  Due  to  Nodal  Loads  Shown  in 
Figure  2.9  _  _ 


1  w 

Elements 

Depth,  y  (m) 

Vertical  Stress, 
ay(kPa) 

Horizontal  Stress, 

CJkPa) 

1,2,3 

0.5 

10 

5 

hti 

1.5 

30 

15 
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These  element  stresses  correspond  to  the  stresses  shown  in  the  stress  profiles  in 
Figure  2.7.  Note  that,  due  to  the  nature  of  the  loading,  shear  stresses  o  are  zero. 


Step  3:  Determination  of  excavation  forces  using  Clough  and  Duncan 
procedure.  To  determine  the  stresses  along  the  excavation  boundaries,  the 
stresses  at  the  nodal  points  2,  3,  6,  and  7  are  first  determined  following  the 
interpolation  procedure  discussed  above.  The  interpolation  is  carried  out  using 
Elements  2, 3,  5,  and  6  as  the  interpolation  elements  and  the  stresses  in  the 
centers  of  these  elements.  The  interpolation  is,  therefore,  carried  out  using  the 
region  connected  by  the  dashed  lines  in  Figure  2.9.  For  illustration,  the  calcula¬ 
tions  are  first  carried  out  for  the  horizontal  stresses  along  the  excavation  bound¬ 
ary  between  Elements  2  and  3  (line  connecting  nodes  3  and  7).  Substituting  the 
horizontal  stresses  and  the  coordinates  of  the  centers  of  Elements  2,  3,  6,  and  5  in 
Equation  2.2  gives  the  following  interpolation  equation: 


1  1.5  1.5 

1  2.5  1.5 

1  2.5  0.5 

1  1.5  0.5 


2.25*1  fo, ' 
3.75  L 
1.25  la,  ’ 
0.75 J  [a4 


Solving  for  {a}  =  {a,  a2  a3  aA  }7  in  Equation  2. 14  gives  {a}  =  {20,  0,-10, 

0}.  Substituting  these  values  in  Equation  2.6,  together  with  the  coordinates  of 
nodes  2,  3, 6,  and  7,  gives  the  following  stresses  at  these  nodal  points: 

Ox  (2)1  1*1  1  2  2*1  f  20 ' 

a,  (3)  1  2  2  4  0 

a, (6)  1212  -10  1 

a,(7)J  [l  1  1  lj[o 

Solving  Equation  2.15  gives 

(2)1  f0 

,0,(3)  I  ° 

0,(«)  10 

,0.(7)]  [lO 


Following  the  same  procedures,  the  vertical  stresses  <jy  at  nodes  2,  3,  6,  and  7 
can  be  obtained  as 
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k(2)1 

'o' 

M3) 

0 

k(6) 

20 

,M7l 

20 

(kPa) 


(2.17) 


Assuming  linear  distribution  of  stresses  along  the  sides  of  the  elements 
between  the  nodes,  the  stresses  along  the  excavation  boundaries  corresponding  to 
#  the  interpolated  nodal  stresses  are  plotted  in  Figure  2.10.  The  boundary  stresses 

are  then  converted  to  nodal  forces  following  the  procedures  shown  in  Figure  2.3. 

The  reversed  equivalents  of  these  nodal  forces  are  also  shown  in  Figure  2.10. 
These  are  the  forces  that  will  be  applied  to  simulate  the  excavation  of  Element  2 
following  the  Clough  and  Duncan  procedure.  Further  discussions  of  the  results  of 
®  the  Clough  and  Duncan  procedure  are  given  in  Section  4.2. 


(a) 


(b) 


Figure  2.10. 


Stress  distribution  (a)  along  the  excavation  boundary  from  the  interpolated  nodal  stresses; 
(b)  reversed  equivalent  nodal  forces  (stresses  are  in  kilopascals;  forces  are  in  kilonewtons) 
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3  Review  of  Other  Excavation 
Models 


As  noted  in  Chapter  2,  several  procedures  for  modeling  excavation  processes 
have  been  developed  since  Clough  and  Duncan  presented  their  excavation 
algorithm  in  1969.  In  this  chapter,  the  main  features  of  the  different  numerical 
models  for  simulating  excavation  processes  are  described  and  reviewed. 


3.1  Summary  of  Models 

The  method  of  using  stresses  or  surface  traction  along  the  excavation  bound¬ 
aries  appears  to  have  been  proposed  first  by  Brown  and  King  (1966).  The  Clough 
and  Duncan  model  itself  was  an  improvement  of  the  models  proposed  by 
Dunlop,  Duncan,  and  Seed  (1968)  and  Chang  (1969).  The  difficulties  associated 
with  the  Clough  and  Duncan  procedure  had  not  been  rigorously  analyzed  until 
Christian  and  Wong  performed  a  detailed  analysis  of  the  different  excavation 
procedures  available  in  1973.  At  that  time,  Clough  and  Duncan  (1969)  and 
Ishihara  (1970)  had  already  demonstrated  that,  for  linear  elastic  materials,  the 
results  of  the  excavation  must  be  independent  of  the  number  of  steps  used  in  the 
excavation. 

Clough  and  Duncan  (1969)  performed  a  numerical  analysis  of  the  effects  of 
number  of  layers  in  simulation  of  one-dimensional  (1-D)  excavation  and  fill 
placement.  Curves  of  error  (in  percent)  in  vertical  displacement  were  developed 
where  the  error  was  defined  as  the  ratio  of  the  vertical  displacement  for  a  given 
number  of  increments  to  the  vertical  displacement  obtained  for  the  solution  using 
layers  of  infinitesimal  thickness.  The  analysis  was  carried  out  for  both  linear  and 
nonlinear  material  behavior.  For  nonlinear  material  response,  the  Young’s  modu¬ 
lus  E  was  modeled  as  a  power  function  of  the  confining  stress  cr3 .  Clough  and 

Duncan  s  results  are  shown  Figure  3.1  in  terms  of  the  values  of  errors  in  simula¬ 
tion  of  small  layers  of  excavation  for  different  values  of  the  exponent  n  used  to 
model  the  dependency  of  E  to  o>3 .  As  can  be  seen,  there  are  no  simulation  errors 

for  linear  elastic  material  with  n  =  0,  which  shows  that  1-D  excavation  is  inde¬ 
pendent  of  the  number  of  layers  used  in  the  simulation. 
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Number  of  Layers 


Figure  3.1.  Values  of  errors  in  simulation  of  small  layer  one-dimensional 
excavation  for  different  values  of  n  (Clough  and  Duncan  1 969) 

Ishihara  (1970)  presented  a  more  formal  and  mathematical  proof  that,  for 
elastic  materials,  the  result  must  be  independent  of  the  number  of  stages  used  in 
general  2-  and  3-D  excavation  problems.  This  follows  from  the  principle  that  the 
response  of  elastic  materials  is  independent  of  the  load  or  displacement  path 
taken  by  the  material.  Uniqueness  of  the  solution  for  excavation  in  elastic 
materials  should,  therefore,  be  one  of  the  main  criteria  in  determining  the  validity 
of  a  numerical  procedure  for  simulating  excavation. 

Christian  and  Wong  (1973)  studied  the  performance  of  three  numerical 
procedures:  (a)  the  Clough  and  Duncan  procedure  of  interpolating  stresses  from 
four  elements  near  the  excavation  boundary;  (b)  the  Dunlop,  Duncan,  and  Seed 
(1968)  procedure  of  averaging  the  stresses  of  the  two  elements  adjacent  to  the 
excavation  boundaiy  ;  and  (c)  using  the  stresses  of  the  excavated  element  as  the 
boundary  stresses.  It  was  found  that  the  results  of  these  three  finite  element 
procedures  could  be  in  serious  error.  In  particular,  the  different  procedures 
yielded  displacements  at  the  top  of  a  vertical  cut  which  were  dependent  on  the 
number  of  steps  of  excavation  in  an  elastic  soil.  For  the  problem  Christian  and 
Wong  (1973)  analyzed,  the  calculated  vertical  displacements  at  the  top  of  the 
excavation  varied  from  -0.6  to  0.8  ft  (-0.2  to  0.2  m)  when  the  number  of 
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excavation  lifts  was  changed  from  1  to  8.  Clearly,  this  discrepancy  shows  the 
difficulty  of  using  the  three  procedures  in  predicting  displacements  from  excava¬ 
tions,  as  the  results  are  heavily  influenced  by  the  number  of  steps  used  in  the 
simulation. 

Christian  and  Wong  (1973)  argued  that  the  errors  are  largely  due  to  the 
inability  of  the  finite  elements  to  model  adequately  the  stress  gradients  at  the  toe 
of  the  excavation.  They  proposed  a  modified  model  and  recommended  that 
simulations  of  excavations  should  be  done  in  as  few  steps  as  possible.  In  their 
modified  procedure,  the  stresses  at  the  excavation  boundaiy  are  extrapolated 
using  a  fifth-  or  sixth-order  polynomial  fitted  through  the  stresses  at  the  element 
centers  in  the  same  horizontal  row  of  elements.  Although  better  results  were 
obtained  with  this  modified  procedure,  there  were  still  differences  in  the  calcu¬ 
lated  displacements  for  different  numbers  of  excavation  steps  used. 

Following  the  work  of  Christian  and  Wong  (1973),  Chandrasekaran  and  King 
(1974)  developed  a  procedure  that  does  not  rely  on  interpolation  of  the  stresses 
along  the  excavation  boundary.  The  procedure  is  based  on  the  individual  calcula¬ 
tions  of  the  excavation  loads  for  each  stage  of  excavation  based  on  earth  pressure 
assumptions  and  the  displacement  functions  used  for  the  elements.  The  indi¬ 
vidual  sets  of  loads  are  calculated  before  the  different  stages  of  excavation  are 
carried  out.  Each  set  of  loads  is  modified  to  account  for  the  changes  in  load  from 
the  previous  excavation  stage.  These  changes  are  equal  to  the  product  of  the  stiff¬ 
ness  matrix  of  the  unexcavated  elements  and  the  displacements  from  the  previous 
stage.  The  procedure  is  shown  to  yield  numerically  identical  results  for  different 
numbers  of  excavation  stages.  Another  attempt  at  solving  the  nonuniqueness 
problem  shown  by  Christian  and  Wong  (1973)  is  the  procedure  developed  by 
Desai  and  Sargand  (1984).  A  hybrid  finite  element  procedure,  in  which  the  finite 
element  degrees-of-freedom  include  both  displacements  and  stresses,  was 
proposed. 

Over  the  years,  a  numerical  solution  that  has  found  wide  acceptance  among 
modelers  is  the  procedure  based  on  calculating  excavation  forces  from  element 
stresses  at  the  Gaussian  points  (or  at  the  element  center  for  constant  stress  ele¬ 
ments).  Clough  and  Mana  (1976)  appeared  to  be  first  to  propose  this  procedure. 
However,  the  most  rigorous  proof  of  the  uniqueness  of  the  solution  based  on  the 
use  of  internal  element  stresses  for  elastic  materials  was  given  by  Ghaboussi  and 
Pecknold  (1984).  A  generalized  formulation  was  presented  in  which  the  excava¬ 
tion  forces  are  treated  in  the  context  of  the  Newton-Raphson  family  of  solution 
techniques  widely  used  in  nonlinear  analysis.  In  the  solution  technique,  unbal¬ 
anced  forces  from  excavation  are  simply  treated  as  part  of  the  force  residuals  that 
need  to  be  balanced  for  a  given  load  or  excavation  stage.  In  addition,  Ghaboussi 
and  Pecknold  accounted  for  the  stresses  from  the  material  self-weight.  The  self¬ 
weight  stresses  were  not  included  in  the  Clough  and  Mana  (1976)  procedure. 

The  validity  of  the  Ghaboussi  and  Pecknold  (1984)  procedure  was  also 
demonstrated  by  Chow  (1985).  Using  this  procedure,  Chow  showed  that  the 
Desai  and  Sargand  (1984)  method  does  not  account  for  the  body  forces.  Brown 
and  Booker  (1986)  performed  further  studies  showing  the  uniqueness  of  the 
solution  using  the  Ghaboussi  and  Pecknold  procedure  for  linearly  elastic 
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materials.  It  was  also  shown  that  the  procedure  does  not  require  a  large  degree  of 
mesh  refinement  to  obtain  accurate  results  for  linear  elastic  materials. 

Boija,  Lee,  and  Seed  (1989)  extended  the  Ghaboussi  and  Pecknold  (1984) 
procedure  to  nonlinear  materials  using  a  variational  formulation  that  accounts  for 
time- varying  problem  domain  and  boundaries.  Independence  of  the  solution  from 
the  number  of  construction  stages  was  shown  not  only  for  elastic  materials  but 
also  for  a  certain  class  of  elastoplastic  materials  with  a  shrinking  elastoplastic 
zone.  Similar  to  the  Ghaboussi  and  Pecknold  procedure,  changes  due  to  excava¬ 
tion  are  treated  as  part  of  the  nonlinear  problem,  although  a  different  numerical 
strategy  is  used  to  resolve  the  force  residuals  from  the  unbalanced  forces. 

Based  on  this  review,  the  methods  can  be  grouped  into  four  main  types  of 
numerical  procedures  for  modeling  excavation  processes,  as  summarized  in 
Table  3.1  and  Figure  3.2.  Chapter  4  presents  an  evaluation  of  these  procedures, 
their  advantages  and  disadvantages,  and  the  requirements  for  implementing  the 
recommended  procedure  in  SOILSTRUCT-ALPHA. 


Table  3.1  1 

Main  Types  of  Numerical  Procedures  for  Simulating  Excavation 
Processes  9 

Numerical  Procedure 

Summary  of  the  Procedure  and  Variants 

Excavation  forces  from 
stresses  along  the 
excavation  boundary 
(Figure  3.2a) 

Originally  developed  by  Brown  and  King  (1966). 

Excavation  boundary  stresses  are  determined  from  the  boundary 
stresses  of  the  excavation  element  (Christian  and  Wong  1973). 
Excavation  boundary  stresses  are  determined  from  the  average  of 
the  element  stresses  of  two  elements  adjacent  to  the  excavation 
boundary  (Dunlop,  Duncan,  and  Seed  1968). 

Stresses  at  the  nodes  of  the  excavated  element  are  determined  from 
interpolation  of  element  stresses  from  four  elements  adjacent  to  the 
excavation  boundary.  Nodal  stresses  are  then  converted  to 
excavation  forces  along  the  excavation  boundary  (Clough  and 

Duncan  1969). 

Excavation  boundary  stresses  are  extrapolated  using  a  fifth-  or  sixth- 
order  polynomial  fitted  through  the  stresses  at  the  element  centers  in 
the  same  horizontal  row  of  elements  (Christian  and  Wong  1973). 

Accumulation  of 
excavation  forces 
(Figure  3.2b) 

Individual  sets  of  excavation  loads  are  calculated  before  the  different 
stages  of  excavation  are  carried  out.  Each  set  of  loads  is  modified  to 
account  for  the  changes  in  the  load  from  the  previous  excavation 
stage.  These  changes  are  equal  to  the  product  of  the  stiffness  matrix 
of  the  unexcavated  elements  and  the  displacements  from  the 
previous  stage  (Chandrasekaran  and  King  1974). 

Hybrid  finite  element 
procedure  (Figure  3.2c) 

Uses  finite  elements  with  displacement  and  stress  degrees-of- 
freedom  (Desai  and  Sargand  1984). 

Excavation  forces  from 
force  residuals 
(Figure  3.2d) 

Excavation  forces  are  calculated  from  stresses  at  the  Gaussian 
points  (Clough  and  Mana  1976). 

Body  forces  are  included  along  with  the  internal  stresses  in  the 
calculation  of  excavation  forces  (Ghaboussi  and  Pecknold  1984, 

Chow  1985,  Brown  and  Booker  1986). 

The  procedure  is  a  generalization  of  the  Newton-Raphson  procedure 
to  minimize  force  residuals  from  unbalanced  (difference  between 
internal  and  external)  forces  (Ghaboussi  and  Pecknold  1984). 

The  procedure  is  generalized  for  nonlinear  materials  behavior  (Borja, 
Lee,  and  Seed  1989). 
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Excavation 
forces  from 
interpolated 
stresses 


Interpolated  stresses 
between  elements 


{Fb}  =  forces  exerted  above  BBB 
on  remaining  soil 
{AFfli}  =  changes  in  {FB}  due  to 
die  first-stage  excavation 


a.  Excavation  forces  from  interpolated 
stresses  along  the  excavation  boundary 


b.  Excavation  forces  accumulated  from 
excavation  stages 


Hybrid  finite  elements 
with  stresses  and 
displacements  at  the 
nodes 


R-Fi \nf-Fp 
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c.  Excavation  forces  from  nodal  stresses 
and  displacement  of  hybrid  finite  elements 


d.  Excavation  forces  calculated  from  force 
residuals  between  internal  and  external 
forces 
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3.2  Methods  Based  on  Force  Residuals 


The  method  of  calculating  excavation  forces  from  the  differences  between 
the  internal  element  forces  and  external  loads  appears  to  be  one  of  the  most 
widely  used  and  perhaps  most  successful  procedure  for  modeling  excavation.  For 
this  reason,  a  detailed  review  of  this  procedure  is  given  in  this  section.  The 
method  is  based  on  the  minimization  of  the  residual  forces  {R}  calculated  from 
the  differences  between  the  internal  forces  {F^}  and  the  external  forces  {Fex,} 
calculated  at  the  nodes  in  a  finite  element  discretization  of  a  stress  problem. 

<31> 


In  a  finite  element  discretization,  {Fmt}  and  {Fext}  are  calculated  from  the 

element  stresses  {c},  body  loads  from  self-weights  {Z>}  =  {0  y}r  where  y  is  the 
unit  weight  (force/length3),  and  surface  stresses  or  tractions  {cr5}  at  the 
boundaries: 


W=t*fw  dr 

V 

(3.2) 

{^-}=  {*}  dV  +  j[N]r{<f}  dS 

(3.3) 

V  s 


where  V  and  S  are,  respectively,  the  volume  and  surface  area  of  the  element,  [Af] 
is  the  global  shape  function  matrix,  and  [5]  is  the  global  strain-displacement 
transformation  matrix. 

For  2-D  stress  conditions,  [5]  is  equal  to 


[B]: 


BN, 

Bx 

0 


0 

BN, 

By 

BN. 

Bx 


(3.4) 


where  /  is  the  node  number  of  the  element. 


So,  in  2-D  plane  stress/plane  strain  problems,  the  internal  forces  are  equal  to 


- 

V 

=J1 


f  BN ,  BN,\  , 

^  ox  oy  j 
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(3.5) 
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Calculations  of  the  internal  and  external  forces  at  element  level  are  illustrated  in 
the  case  of  the  isoparametric  four-noded  quadrilateral  element  as  shown  in 
Appendix  A. 

The  internal  and  external  loads  are  summed  for  all  elements  connected  to  a 
node  from  the  global  internal  and  external  force  vectors  {Fmt}  and  {F^,}  . 

During  excavation,  both  of  these  force  vectors  will  change,  and  the  differences 
between  the  two  force  vectors  correspond  to  the  excavation  loads  at  the  excava¬ 
tion  boundaries.  The  changes  in  the  internal  forces  are  due  to  the  removal  of  the 
element  stresses  of  the  excavated  elements.  For  the  deleted  elements,  the  internal 
stresses  {a}  in  Equation  3.2  are  set  to  zero;  thus,  deleted  elements  no  longer 
contribute  to  the  internal  forces. 

The  changes  in  the  external  forces  are  due  to  the  removal  of  the  nodal  forces 
for  the  nodes  within  the  excavation  and  to  the  reduction  of  the  nodal  forces  in  the 
nodes  along  the  excavation  boundaries.  For  the  deleted  elements,  the  body  loads 

{6}  and  the  surface  tractions  on  the  boundary  of  the  deleted  elements  {cr,J  in 

Equation  3.3  are  set  to  zero;  thus,  deleted  elements  no  longer  contribute  to  the 
external  forces. 

The  difference  {F„,}  -{F’nt}  will  cause  deformations  {Am}  ,  which  are 
computed  from  the  global  force-displacement  relation: 

[F]{Am}={F„,  }-{/;„,}  (3.6) 

where  [F]  =  J[fi]7  [Z)][Z?]c/F  is  the  global  stiffness  matrix  ([/)]  is  referred  to  as 

V 

either  the  constitutive  or  stress-strain  matrix).  The  deformations  {Am}  are  the 
deformations  induced  by  the  excavation. 

One  approach  to  apply  the  method  of  force  residuals  is  illustrated  in  a 
flowchart  in  Figure  3.3. 

Consistent  nodal  point  forces  due  to  the  self-weight  of  each  element  are 
computed  in  Equation  3.3  using  the  following  {&} -matrix: 

W  =  {°}  (3-7) 

where  y  is  the  unit  weight  of  the  soil.  If  the  element  is  rectangular,  the  volume 

integral  of  [A]  {0  y}  in  Equation  3.3  predicts y-direction  loads  at  each  node, 

each  equal  to  one  quarter  of  the  total  forces  (i.e.,  the  total  weight)  on  the  element. 
This  is  regardless  of  the  orientation  of  the  element  with  respect  to  the  y-axis 
(Cook,  Malkus,  and  Plesha  2002,  p  114). 


Chapter  3  Review  of  Other  Excavation  Models 


Perform  gravity  turn-on  analysis  to 
determine  the  initial  in  situ  stresses. 
Set  displacements  to  zero. 


Figure  3.3.  Flowchart  for  excavation  modeling  using  the  method  of  force 
residuals 
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3.3  Example  Calculation 


The  problem  analyzed  using  the  Clough  and  Duncan  procedure  in 
Section  2.2  will  be  reanalyzed  to  illustrate  the  use  of  the  procedure  described 
above.  For  convenience,  the  finite  element  mesh  (Figure  2.9)  is  reproduced 
below  as  Figure  3.4. 


V 


Figure  3.4.  Nodal  loads  from  self-weight  (forces  shown  are  in  kilonewtons) 


Steps  1  and  2:  Gravity  turn-on  analysis  and  results  of  analysis 

These  steps  are  the  same  as  those  used  in  the  Clough  and  Duncan  procedure. 
At  the  end  of  this  step,  the  nodal  load  vectors  due  to  element  self-weights  are  as 
shown  in  Figure  3.4. 

Step  2:  Calculation  of  equivalent  internal  element  forces  from 
internal  stresses 

These  are  the  nodal  forces  that  are  equivalent  to  the  internal  stresses.  The 
authors  of  this  report  will  refer  to  them  as  the  internal  element  forces.  This  step 
corresponds  to  the  calculation  of  the  equivalent  element  nodal  forces  ( Fxi  )  and 

(fyt ).m for  each  element  using  Equation  3.5  corresponding  to  the  stresses  ox,  oy, 
and  cxy  in  the  element.  Recall  that  the  element  stresses  due  to  the  self-weight 
loads  in  Elements  1, 2,  and  3  are  oy=  10  kPa  and  ox  =  5  kPa,  and  the  corre¬ 
sponding  stresses  for  Elements  4,  5,  and  6  are  oy=  30  kPa  and  <^=15  kPa.  For 
all  elements,  Oxy  =  0.  The  calculations  of  element  nodal  forces  are  carried  out  for 
all  the  element  nodes  /  (=  1  to  4  for  the  four-noded  constant  stress  element). 
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In  the  case  of  constant  stress  elements  this  operation  is  trivial,  as  the  element 
forces  are  simply  based  on  force  equilibrium.  A  constant  stress  distribution  will 
yield  an  equal  distribution  of  nodal  forces  among  the  nodes.  This  is  shown  in 
Figure  3.5  below  in  the  calculation  of  the  horizontal  nodal  forces  at  the  two 
nodes  on  the  right  side  of  Element  2.  The  horizontal  nodal  forces  are  simply  what 
are  needed  to  balance  the  horizontal  stress  a,  =  5  kPa  in  the  elements  based  on 
the  free-body  diagram  shown. 


Figure  3.5.  Calculation  of  reversed  internal  element  nodal  forces  based  on  force 
equilibrium 

The  calculation  of  the  internal  element  forces  using  Equations  3.2  and  3.5  is 
also  illustrated  in  the  Appendix  A.  For  a  general  loading  condition,  the  four- 
noded  isoparametric  element  discussed  in  Appendix  A  does  not  in  general  have  a 
constant  strain  nor  constant  stresses  throughout  the  element.  It  can  be  shown  that 
for  the  loading  condition  (i.e.,  gravity  turn-on)  applied  to  the  Figure  3.4  (homog¬ 
enous,  elastic)  regular  mesh,  the  four-noded  isoparametric  element  has  constant 
stress  throughout  each  element.  For  the  isoparametric  element,  Equation  A.22 
results  in  the  equivalent  nodal  point  forces  given  in  Figure  3.5.  Appendix  A 
discusses  the  different  methodologies  for  computing  the  equivalent  internal 
forces  from  the  stress  field  within  the  element.  The  horizontal  forces  in  all  the 
nodes  of  the  element  are  shown  to  be  equal  to  2.5  kN,  as  is  obtained  in  Figure  3.5 
and  also  from  Equation  A.22  with  a  uniform  stress  field.  The  complete  element 
forces  corresponding  to  the  element  stresses  of  oy  =  10  kPa  and  ax=  5  kPa  in 

Elements  1,  2,  and  3,  and  oy=  30  kPa  and  ax=  15  kPa  in  Elements  4,  5,  and  6 

are  shown  in  Figure  3.6. 

Step  3:  Calculation  of  internal  nodal  forces  from  equivalent 
internal  element  forces 

In  this  step,  the  internal  nodal  forces  from  all  the  element  stresses  are 
obtained.  This  is  done  by  summing  the  internal  element  forces  from  all  the 
elements  connected  to  each  node  after  the  excavation.  The  results  from  this  step 
are  shown  in  Figure  3.7.  Note  that  only  the  forces  along  the  excavation  bound¬ 
aries  need  to  be  calculated,  as  shown  in  Figure  3.7.  It  can  be  demonstrated  that 
the  summation  of  internal  and  external  forces  in  the  nodes  that  do  not  coincide 
with  the  excavation  boundaries  will  be  zero,  and  these  nodes  will  be  in  equili¬ 
brium.  The  element  internal  forces  are  shown  for  both  the  excavated  element 
(Element  2)  and  the  remaining  soil.  The  internal  forces  for  Element  2  are  the 
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Figure  3.6.  Reversed  internal  nodal  forces  from  element  stresses 


same  as  those  in  Figure  3.6.  For  the  remaining  soil,  the  nodal  forces  are  calcu¬ 
lated  from  element  force  contributions  from  all  elements  connected  to  a  node. 

For  example,  consider  the  upward  vertical  force  of  25  kN  for  node  6  in  Fig¬ 
ure  3.7b.  The  elements  connected  to  node  6  after  the  excavation  are  1, 4,  and  5. 
The  upward  vertical  force  of  25  kN  for  node  6  is  obtained  from  summation  of  the 
downward  5  kN  from  Element  1,  the  upward  force  of  15  kN  for  Element  4,  and 
another  upward  force  of  1 5  kN  for  Element  5.  Similar  operations  are  carried  out 
for  all  the  nodes  for  both  horizontal  and  vertical  forces. 


Figure  3.7.  Reversed  total  internal  element  forces  along  the  excavation  boundaries  (forces  are  in 
kilonewtons) 
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Step  4:  Calculation  of  external  loads 

The  external  nodal  loads  due  to  the  self-weight  on  the  excavated  element  and 
the  remaining  soil  after  the  excavation  are  calculated  similarly  to  the  loads  in 
Figure  3.4.  The  new  external  nodal  loads  after  excavation  are  shown  in  Fig¬ 
ure  3.8.  For  the  excavated  Element  2,  the  external  loads  are  simply  equal  to 
one  quarter  of  the  weight  of  the  element  (20  kN/4  =  5  kN)  equally  distributed 
among  the  four  nodes.  For  the  remaining  soil,  the  main  difference  from  Fig¬ 
ure  3.4  is  that  Element  2  no  longer  contributes  to  the  external  loads.  Thus,  the 
external  loads  at  nodes  2,  3,  6,  and  7  will  have  to  be  modified  to  account  for  the 
deletion  of  Element  2.  The  external  loads  in  the  other  nodes  not  in  the  excavation 
boundary  will  be  unchanged.  For  example,  the  nodal  force  in  node  2,  which  was 
10  kN  in  Figure  3.4,  is  now  only  5  kN.  This  force  comes  from  one  quarter  of  the 
weight  of  Element  1,  the  only  element  left  connected  to  node  2  after  the 
excavation. 


Figure  3.8.  Total  external  nodal  loads  along  the  excavation  boundaries  (forces  are  in  kilonewtons) 


Step  5:  Calculate  unbalanced  excavation  forces 

After  the  gravity  turn-on  analysis,  the  internal  forces  and  external  forces  are 
equal  and  the  system  is  in  equilibrium.  However,  the  excavation  of  Element  2 
will  create  unbalanced  forces  from  the  differences  in  the  new  external  and 
internal  forces.  Adding  the  reversed  internal  forces  (Figure  3.7)  to  the  corre¬ 
sponding  external  loads  from  the  self-weights  (Figure  3.8)  gives  the  net  unbal¬ 
anced  excavation  forces  along  the  excavation  boundaries  (Figure  3.9).  These  are 
the  forces  that  will  cause  displacements  due  to  excavation. 

Note  the  correspondence  between  the  excavation  forces  on  the  excavated 
element  and  the  remaining  soil.  The  forces  are  of  equal  magnitudes  but  are 
opposite  in  sign.  This  means  that  excavation  forces  can  be  calculated  using  either 
the  excavated  element(s)  or  the  remaining  soil,  and  the  results  will  be  the  same. 
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Figure  3.9.  Excavation  forces  along  the  excavation  boundaries  (forces  are  in  kilonewtons) 

Note  also  the  differences  between  the  horizontal  excavation  forces  shown  in 
Figure  3.9b  and  those  obtained  from  the  Clough  and  Duncan  procedure  in 
Figure  2.10.  The  results  of  the  procedure  based  on  minimization  of  residual 
forces  are  discussed  further  in  Section  4.2. 
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4  Assessment  of  Numerical 
Models  for  Excavation 


This  chapter  provides  a  detailed  comparison  of  the  numerical  procedures  for 
simulating  excavation  presented  in  Chapter  3.  The  different  numerical  procedures 
are  evaluated  in  terms  of  uniqueness  of  the  solution,  satisfaction  of  force  equi¬ 
librium,  and  convergence  of  the  solution.  Comparisons  are  also  made  in  terms  of 
practical  issues,  such  as  ease  of  implementation  in  SOILSTRUCT-ALPHA. 


4.1  Comparisons  of  Models 

Four  types  of  numerical  models  for  simulating  excavation  were  identified 
from  the  review  performed  in  Chapter  3,  and  the  main  features  of  these  numerical 
models  were  summarized  in  Table  3.1.  Based  on  this  review,  the  advantages  and 
disadvantages  of  the  four  methods  were  also  identified.  A  brief  summary  of  the 
major  advantages  and  disadvantages  of  the  methods  is  given  in  Table  4.1. 


Table  4.1 

Advantages  and  Disadvantages  of  Different  Types  of  Numerical 
Procedures  for  Simulating  Excavation  Processes 

Numerical 

Procedure 

Advantages/Disadvantages 

Excavation  forces 
from  stresses  along 
the  excavation 
boundary 

Attempts  to  account  for  stress  concentrations  at  sharp  excavation  corners 
and  between  elements  of  different  sizes. 

Does  not  guarantee  force  equilibrium  between  internal  element  forces  and 
external  loads. 

Results  are  dependent  on  the  excavation  sequence,  and  do  not  guarantee 
a  unique  solution  for  elastic  problems. 

Accumulation  of 
excavation  forces 

Solution  is  independent  of  the  number  of  excavation  stages  for  elastic 
problems. 

Difficult  to  program,  particularly  for  complicated  excavation  sequences  and 
geometries. 

Requires  large  computer  memory  to  store  loads  from  previous  excavation 
sequences. 

Hybrid  finite  element 
procedure 

Solution  is  independent  of  the  number  of  excavation  stages  for  elastic 
problems. 

Uses  nonstandard  element  formulations. 

Difficult  to  program  and  requires  more  computer  resources. 

Excavation  forces 
from  force  residuals 

Consistent  formulation,  guaranteeing  force  equilibrium  between  internal 
element  forces  and  external  loads. 

Results  are  independent  of  the  excavation  sequence  for  elastic  materials. 
Uses  standard  finite  element  calculations  (e.g.,  use  of  the  [B]-matrix  and 
Gaussian  quadrature  to  calculate  internal  elements  forces). 

Can  be  extended  to  nonlinear  problems  in  a  consistent  manner. 
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As  mentioned  in  Chapter  3,  although  the  Clough  and  Duncan  (1969) 
approach  was  one  of  the  earliest  procedures  to  be  widely  used  to  analyze  exca¬ 
vation  problems,  Christian  and  Wong  (1973)  later  recognized  potential  problems 
with  this  approach.  One  difficulty  of  using  the  Clough  and  Duncan  procedure  is 
that  the  results  of  the  modeling  are  strongly  dependent  on  the  number  of  steps 
used  to  simulate  the  excavation.  The  method  they  suggested  to  minimize  the 
dependency  of  results  on  the  number  of  excavation  steps  was  to  model  the  exca¬ 
vation  in  as  few  steps  as  possible.  However,  the  suggested  method  goes  against 
the  philosophy  in  finite  element  analysis  of  achieving  a  converging  solution  as 
the  problem  domain  is  refined  (both  in  size  and  the  load  increment).  It  also  poses 
a  difficulty  for  nonlinear  problems  where  the  solution  calls  for  the  application  of 
small  load  increments  at  a  time. 

Christian  and  Wong  (1973)  believed  that  the  discrepancies  in  the  results  were 
due  to  the  inability  of  finite  elements  to  model  high  stress  concentrations  at  the 
toe  of  the  excavation.  However,  even  with  an  improved  method  to  extrapolate 
stresses  along  the  excavation  boundaries  using  a  fifth-order  polynomial  (in  com¬ 
parison  to  the  bilinear  interpolation  used  in  the  Clough  and  Duncan  procedure), 
the  results  were  still  dependent  on  the  number  of  excavation  steps  used.  In  fact, 
the  difficulties  associated  with  the  Clough  and  Duncan  procedure  can  be  attribu¬ 
ted  to  the  attempt  at  arriving  at  an  accurate  representation  of  the  stress  concen¬ 
trations  caused  by  excavation.  It  was  argued,  quite  reasonably,  that  an  accurate 
representation  of  the  stresses  along  the  excavation  boundaries  is  essential  in 
obtaining  an  accurate  solution  of  the  excavation  process.  Unfortunately,  their 
method  can  produce  boundary  stresses  that  are  not  in  equilibrium  with  the  ele¬ 
ment  stresses  and  self-weights  of  the  excavated  elements. 

Although  three  other  methods  of  excavation  modeling  were  identified,  two  of 
these  methods  have  fallen  out  of  use  in  recent  years.  The  method  of  Chandrase- 
karan  and  King  (1974)  is  based  on  the  accumulation  of  loads  from  previous 
stages  of  excavation  and  requires  the  excavated  region  to  be  divided  into  differ¬ 
ent  stages.  The  loads  at  the  nodes  corresponding  to  the  boundaries  dividing  the 
excavation  stages  must  be  recalculated  and  stored  at  every  stage  of  the  excava¬ 
tion.  This  method  is  cumbersome,  difficult  to  program  (particularly  for  compli¬ 
cated  excavation  sequences  and  geometries),  and  uses  a  lot  of  computer  memory. 

Another  method  that  has  not  found  wide  application  is  that  developed  by 
Desai  and  Sargand  (1984).  Again,  this  method  is  an  outcome  of  the  belief  that 
stresses  along  the  excavation  boundary  must  be  accurately  represented  to  model 
excavation  accurately.  To  accomplish  this,  Desai  and  Sargand  used  a  hybrid 
finite  element  formulation  where  stresses,  in  addition  to  the  usual  displacements, 
are  treated  as  nodal  unknowns.  The  result  of  this  formulation,  however,  is  an 
increase  in  complication  of  the  finite  element  solution.  Hybrid  finite  element 
formulations  are  more  difficult  to  program  and  require  more  computational 
resources  due  to  the  increased  number  of  nodal  unknowns.  Moreover,  it  is  not 
clear  how  the  Desai  and  Sargand  procedure  can  be  extended  to  account  for  body 
loads  from  element  self-weights. 

Of  the  four  methods  reviewed,  it  appears  that  the  most  rigorous  and  analyti¬ 
cally  robust  method  is  the  one  based  on  balancing  residual  forces  from  the  differ¬ 
ences  between  the  internal  and  external  loads  in  a  loaded  body.  This  method  is 
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now  the  most  widely  used  and  accepted  procedure  for  excavation  modeling  (e.g., 
Ghaboussi  and  Pecknold  1984;  Borja,  Lee,  and  Seed  1989;  Langer  and 
Stockmann  1985;  Morrison  and  Duncan  1996).  The  objective  of  this  method  is 
first  and  foremost  to  balance  any  unbalanced  loads  created  by  the  excavation. 

The  method  is  based  on  the  recognition  that,  prior  to  excavation,  a  body  of  soil  is 
in  a  state  of  equilibrium  with  the  internal  stresses  (stresses  in  the  elements)  just 
balancing  the  external  loads  (from  body  weights  and  external  surface  tractions). 
Excavation  disturbs  this  equilibrium  when  stressed  elements  are  taken  out 
(changing  the  internal  loads)  and  the  loads  from  the  self-weights  are  reduced  by 
the  excavation  (changing  the  external  loads).  By  ensuring  that  force  equilibrium 
is  satisfied  at  every  stage  of  calculation,  the  method  based  on  minimization  of 
residual  forces  guarantees  a  unique  solution  for  any  number  of  excavation  stages 
for  an  elastic  material.  This  was  proven  formally  by  Ghaboussi  and  Pecknold 
(1984). 

Aside  from  rigorously  satisfying  force  equilibrium,  the  method  based  on 
force  residuals  has  other  distinct  advantages  over  the  other  methods.  One  major 
advantage  is  that  the  method  requires  quantities  that  can  be  directly  obtained 
from  standard  finite  element  calculations.  This  is  particularly  true  for  the  calcula¬ 
tion  of  the  forces  from  internal  stresses  based  on  Equation  3.2.  Calculation  of 
element  stresses  is  routinely  done  in  all  finite  element  codes.  Also,  the  [  B\  - 

matrix  is  already  available  from  the  calculation  of  the  stiffness  matrix,  and  it  can 
be  stored  for  later  use  or  recalculated  when  needed  in  the  evaluation  of  the 
internal  nodal  forces.  Similarly  for  the  recalculation  of  the  external  loads,  the 
new  external  loads  may  either  be  prescribed  manually  or  calculated  from  Equa¬ 
tion  3.3.  The  numerical  integration  required  in  Equations  3.2  and  3.5  can  be  done 
using  routines  for  generating  element  matrices.  In  fact,  most  of  the  calculations 
are  standard  and  can  be  done  at  the  element  level,  and  later  assembled  following 
standard  finite  element  assembly  procedures. 

Another  major  advantage  is  that  the  method  based  on  balancing  of  force 
residuals  lends  itself  naturally  to  nonlinear  problems,  as  shown  by  Ghaboussi  and 
Pecknold  (1984)  and  Borja,  Lee,  and  Seed  (1989).  In  the  case  of  nonlinear  prob¬ 
lems,  the  force  residuals  are  due  to  the  loads  from  differences  in  stresses  between 
the  correct  and  the  current  solution  (the  so-called  initial  stresses  in  nonlinear 
analysis).  In  fact,  in  the  formulations  of  Ghaboussi  and  Pecknold  (1984)  and 
Borja,  Lee,  and  Seed  (1989),  changes  in  the  problem  geometry  from  excavation 
and  addition  of  soil  mass  (e.g.,  from  embankments)  are  simply  treated  as  part  of 
the  nonlinear  problem.  In  this  manner,  excavation  modeling  can  take  advantage 
of  available  methods  to  solve  nonlinear  problems. 


4.2  Performance  of  Numerical  Procedures  for 
Excavation  Modeling 

Based  on  the  above  comparisons,  it  is  apparent  that  the  most  viable  and 
realistic  alternative  to  the  procedure  based  on  interpolation  of  stresses  along  the 
excavation  boundary,  in  particular  the  Clough  and  Duncan  procedure  used  in 
SOILSTRUCT-ALPHA,  is  the  one  based  on  the  minimization  of  force  residuals. 
The  other  methods  based  on  accumulation  of  excavation  loads  and  the  use  of 
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hybrid  finite  element  formulation  have  not  gained  much  acceptance,  and  are 
computationally  more  difficult  to  implement  and  use.  As  mentioned,  the  method 
based  on  achieving  force  equilibrium  between  internal  and  external  forces  is  cur¬ 
rently  the  most  accepted  and  widely  used  method  for  modeling  excavation.  In  the 
following  sections,  further  comparisons  are  made  between  the  procedure  based 
on  interpolation  of  boundary  stresses  and  the  method  based  on  force  residuals.  As 
the  Clough  and  Duncan  procedure  is  used  to  interpolate  the  stresses  along  the 
excavation  boundary  in  SOILSTRUCT,  the  comparisons  here  will  be  made  using 
the  Clough  and  Duncan  procedure.  The  comparisons  and  performance  analyses 
are  made  to  show  whether  and  how  the  two  methods  meet  three  fundamental 
requirements:  (a)  the  solution  must  be  unique  for  path-independent  linear  prob¬ 
lems,  (b)  force  equilibrium  must  be  satisfied  at  every  stage  of  the  calculation,  and 
(c)  the  solution  must  be  convergent  with  refinement  of  the  finite  element 
discretization. 


4.2.1  Uniqueness  of  solution  for  homogenous  elastic  regime 

As  discussed  in  Chapter  3,  the  uniqueness  of  solution  for  excavation  in 
homogenous  elastic  media  has  been  demonstrated  by  Clough  and  Duncan  (1969) 
for  one-dimensional  (1-D)  excavation  and  fill  placement,  and  more  formally  by 
Ishihara  (1970)  for  general  2-  and  3-D  excavation.  The  simple  problem  shown  in 
Figure  4.1a  involving  excavation  in  an  elastic  soil  will  be  used  to  illustrate 
whether  the  methods  based  on  force  residuals  and  the  Clough  and  Duncan  pro¬ 
cedure  yield  unique  solutions  for  elastic  materials  independent  of  the  number  of 
excavation  steps.  The  problem  consists  of  16  four-noded  isoparametric  quadri¬ 
lateral  finite  elements  with  the  same  elastic  properties  for  all  elements.  The  soil 
has  a  unit  weight  of  120  pcf,  a  Young’s  modulus  of  E  =  l-lO6  psf  and  Poisson’s 
ratio  ofv  =  0.3.  Each  element  is  10  by  10  ft  (3  by  3  m)  in  size.  The  objective  is  to 
simulate  the  displacements  of  the  soil  after  Elements  1, 2,  5,  and  6  are  deleted 
from  the  problem  domain.  Roller  boundaries  are  used  on  the  sides  and  the  base  of 
the  model  to  prevent  displacements  normal  to  the  boundaries  but  allow  displace¬ 
ments  parallel  to  the  boundaries.  The  first  step  in  the  modeling  is  the  gravity 
turn-on  analysis  to  determine  the  initial  stresses  from  the  element  self-weights. 
This  is  followed  by  the  deletion  of  the  elements  in  the  excavated  region. 

Uniqueness  of  solution  will  be  investigated  first  in  the  case  of  the  method  of 
force  residuals.  Figure  4.1b  shows  the  deformed  mesh  (with  displacements  exag¬ 
gerated  lOx)  after  the  excavation  is  made  in  one  stage  (i.e..  Elements  1,  2,  5,  and 
6  were  taken  out  at  the  same  time).  The  most  prominent  effect  of  the  excavation 
is  the  heaving  at  the  base  of  the  excavation,  with  the  maximum  upward  move¬ 
ment  of  the  base  equal  to  0.401  ft  (0.12  m).  There  are  only  minor  lateral  move¬ 
ments  along  the  vertical  edge  of  excavation. 

Figures  4.2a  and  b  show  the  results  for  the  case  of  a  two-step  excavation.  The 
nodal  forces  shown  in  these  two  figures  are  discussed  below.  In  the  first  stage  of 
the  excavation,  only  the  first  layer  (consisting  of  Elements  1  and  5)  was  deleted. 
The  deformed  mesh  after  this  first  stage  of  excavation  is  shown  in  Figure  4.2a.  In 
the  second  stage,  Elements  2  and  6  (corresponding  to  the  second  row)  are 
deleted.  The  deformed  mesh  after  the  second  stage  excavation  is  shown  in 
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Figure  4.1 .  Method  of  force  residuals  using  four-nodal  isoparametric  quadrilateral 
elements — single-stage  excavation 
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Figure  4.2b.  The  maximum  heave  at  the  base  of  the  excavation  is  also  0.401  ft. 
Comparison  of  Figures  4.1b  and  4.2b  shows  that  the  two  figures  are  identical, 
demonstrating  that  the  results  of  one-  and  two-stage  excavation  are  identical  for 
this  case. 

To  illustrate  further  the  uniqueness  of  excavation  solutions  for  elastic  mate¬ 
rials  using  the  method  of  force  residuals,  the  same  problem  shown  in  Figure  4.1a 
is  reanalyzed  using  higher  order  eight-noded  isoparametric  quadrilateral  elements 
(Figure  4.3a).  Due  to  the  higher  order  of  the  element  used,  the  mesh  shown  in 
Figure  4.3a  has  twice  as  many  nodes  as  in  Figure  4.1a.  The  same  properties  and 
dimensions  are  used  as  those  in  Figure  4. 1 .  The  results  of  the  one-  and  two-stage 
excavation  are  shown  in  the  deformed  meshes  with  exaggerated  deformations  in 
Figure  4.3b  and  4.4,  respectively.  Again,  identical  results  are  obtained  from  the 
single-  and  the  two-stage  excavations.  The  maximum  heave  at  the  base  of  the 
excavation  is  0.399  ft  for  both  Figures  4.3b  and  4.4b,  which  is  only  slightly 
smaller  than  the  maximum  heave  of  0.401  ft  from  the  analysis  using  four-noded 
elements. 

The  same  problem  shown  in  Figure  4.4  is  now  analyzed  using  the  method 
based  on  the  Clough  and  Duncan  interpolation  of  stresses  along  the  excavation 
boundary.  The  result  of  the  one-stage  simulation  is  shown  in  a  deformed  mesh  in 
Figure  4.5a.  The  maximum  heave  at  the  base  of  the  excavation  is  equal  to 
0.406  ft,  which  is  only  slightly  higher  than  the  corresponding  value  of  0.401  ft 
for  the  method  of  force  residuals.  A  comparison  of  this  result  with  that  obtained 
from  the  method  of  force  residuals  is  shown  in  Figure  4.5b.  As  can  be  seen,  there 
is  very  little  difference  in  the  results  from  the  two  methods,  showing  that  the  one- 
stage  Clough  and  Duncan  procedure  can  yield  comparable  results  to  the  method 
of  force  residuals,  at  least  for  the  simple  problem  analyzed. 

The  results  of  the  two-stage  excavation  using  the  Clough  and  Duncan  pro¬ 
cedure  are  shown  in  Figure  4.6.  The  deformed  mesh  at  the  end  of  the  first  stage 
excavation  is  shown  in  Figure  4.6a,  and  for  the  second  stage  in  Figure  4.6b.  (The 
nodal  forces  shown  in  Figures  4.6a  and  4.6b  are  discussed  below.)  The  maximum 
heave  at  the  base  of  the  excavation  is  about  0.419  ft,  which  is  slightly  larger  than 
that  for  the  single-stage  excavation.  A  comparison  of  the  deformed  mesh  at  the 
end  of  the  complete  excavation  from  the  single-  and  two-stage  excavation  using 
the  Clough  and  Duncan  procedure  is  shown  in  Figure  4.7.  As  can  be  seen,  the 
results  are  different  with  the  two-stage  excavation,  showing  more  heave  at  the 
base  of  the  excavation  and  at  the  top  of  the  vertical  edge  of  the  excavation  than 
the  one-stage  excavation.  This  result  is  counter  to  the  requirement  that  the  results 
be  independent  of  the  number  of  excavation  stages  for  linear  elastic  materials. 

It  is  expected  that  differences  between  multistage  excavation  and  single-stage 
excavation  will  increase  with  increasing  number  of  excavation  stages  in  the 
Clough  and  Duncan  procedure.  This  is  due  to  the  accumulative  nature  of  the  error 
in  the  magnitude  and  distribution  of  excavation  forces  in  their  procedure  (Desai 
and  Sargand  1984).  One  approach  to  reduce  inaccuracies  in  the  Clough  and 
Duncan  procedure  is  to  perform  as  few  excavation  steps  as  possible,  and  prefer¬ 
ably  use  only  a  single  stage  of  excavation.  The  results  given  above  show  that,  at 
least  for  simple  geometries,  the  single-stage  excavation  approach  can  yield  simi¬ 
lar  results  to  those  obtained  from  the  method  of  force  residuals.  Christian  and 
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Figure  4.3.  Method  of  force  residuals  using  eight-nodal  isoparametric 
quadrilateral  elements — single-stage  excavation 
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y  axis  »  y  axis 


Figure  4.5.  Comparison  the  results  of  a  single-stage  excavation  using  the  Clough 
and  Duncan  procedure  and  the  Method  of  Force  Residuals 
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Figure  4.6.  Clough  and  Duncan  procedure — two-stage  excavation  (deformations 
magnified  lOx)  (nodal  forces  shown  are  in  pounds) 
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Figure  4.7.  Comparison  of  the  results  of  a  single-stage  (thin  lines)  and  a  two- 
stage  (thick  lines)  excavation  using  the  Clough  and  Duncan 
procedure  (deformations  magnified  lOx) 


Wong  (1973)  also  provided  numerical  results  to  support  the  use  of  few  excava¬ 
tion  stages.  It  is  should  be  noted,  however,  that  this  recommendation  runs  counter 
to  the  requirement  of  using  small  load  increments  for  problems  with  significant 
material  nonlinearity.  Also,  the  nature  of  the  construction  (e.g.,  excavations  using 
tieback  construction)  can  set  a  limit  on  the  minimum  number  of  excavation 
stages  that  can  be  simulated. 

The  results  of  the  different  simulations  are  summarized  in  Table  4.2  in  terms 
of  the  maximum  heave  at  the  base  of  the  excavation.  The  method  offeree  residu¬ 
als  gave  displacements  of  0.401  ft  for  both  the  one-  and  the  two-stage  excavation 
using  four-noded  elements.  The  corresponding  results  using  the  eight-noded 
elements  are  only  slightly  less  at  0.399  ft,  showing  very  minimal  dependency  of 
the  results  from  the  level  of  discretization.  The  results  of  the  Clough  and  Duncan 
procedure  gave  a  maximum  heave  of  0.406  ft  for  the  single-stage  excavation  and 
0.419  ft  for  the  two-stage  excavation. 


|  Table  4.2 

1 

|  Calculated  Maximum  Heave  at  End  of  Excavation  1 

B  Model 

One-Stage  Excavation 

Two-Stage  Excavation  i 

|  Method  of  force  residuals 

1  4-noded  elements 

3  8-noded  elements 

0.401  ft 

0.399  ft 

— .  11  1  sss=s 

0.401  ft 

0.399  ft 

|  Clough  and  Duncan  procedure 

1  4-noded  elements 

0.406  ft 

0.419  ft  1 

|  Note:  To  convert  feet  to  meters,  multiply  by  0.3048.  | 
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The  larger  heave  at  the  base  of  the  excavation  and  of  the  original  ground 
surface  at  the  top  of  the  vertical  excavation  boundary  for  the  two-stage  excava¬ 
tion  compared  to  the  single-stage  excavation  (Figure  4.7)  is  attributed  to  the 
^  larger  vertical  excavation  forces  in  Figure  4.6b.  This  will  be  discussed  further  in 

Section  4.2.2  in  connection  with  force  equilibrium.  The  results  shown  in  Figure 
4.7  illustrate  that  the  Clough  and  Duncan  procedure  does  not  yield  a  unique  result 
independent  of  the  number  of  excavation  stages,  as  expected  for  linear  elastic 
materials. 

•  The  element  internal  forces  were  calculated  using  stresses  evaluated  at  four 

Gaussian  quadrature  points.  These  are  the  internal  forces  used  in  the  force  resid¬ 
ual  method  shown  in  Figures  4. 1-4.4.  Note  also  that  there  are  horizontal  forces  at 
the  lateral  vertical  boundaries  due  to  the  gravity  turn-on.  However,  these  forces 
are  not  active  since  the  nodes  along  the  lateral  boundaries  are  fixed  along  the  hor¬ 
izontal  direction.  These  forces  are,  therefore,  not  shown  in  Figures  4.2  and  4.6. 


4.2.2  Force  and  moment  equilibrium 

The  lack  of  moment  equilibrium  in  the  Clough  and  Duncan  procedure  for 
constant  stress  elements  can  be  shown  by  recalling  the  example  problem  that  was 
analyzed  in  Section  2.2.  At  the  end  of  the  calculation,  the  stresses  along  the  exca¬ 
vation  boundaries  and  the  equivalent  nodal  forces  (but  opposite  in  sign)  from 
these  stresses  are  as  shown  in  Figure  4.8  below.  It  can  be  seen  that  equilibrium  of 
forces  is  satisfied  along  the  vertical  direction,  as  shown  in  Figure  4.9,  since  the 
nodal  force  of  2(10)  kN  equals  half  the  weight  of  the  element  (10  kN)  plus  the 
effect  of  vertical  stress  at  the  element  center  of  (10  kPa)(l  m2)  =  10  kN. 

In  the  horizontal  direction,  the  horizontal  forces  required  to  satisfy  force  and 
moment  equilibrium  at  the  top  and  bottom  nodes  are  2.5  kN  for  a  constant  stress 
element.  However,  the  nodal  forces  are  equal  to  1.67  kN  at  the  top  nodes  and 
3.33  kN  at  the  bottom  nodes  (in  Figure  4.8b)  according  to  the  Clough  and 
Duncan  procedure.  Horizontal  force  equilibrium  is  achieved  only  on  the  average, 
as  the  average  horizontal  nodal  force  is  equal  to  (1.67  +  3.33)/2  =  2.5  kN.  How¬ 
ever,  moment  equilibrium  will  not  be  satisfied  for  a  constant  stress  element. 

In  comparison,  the  excavation  nodal  forces  obtained  by  minimizing  the 
residual  forces,  which  are  reproduced  in  Figure  4.10,  satisfy  force  and  moment 
equilibrium  both  in  the  horizontal  and  vertical  directions. 

More  convincing  evidence  of  the  lack  of  force  equilibrium  in  the  Clough  and 
Duncan  procedure  can  be  obtained  from  the  results  of  the  finite  element  analysis 
carried  out  in  Section  4.2.1,  particularly  by  comparing  the  excavation  nodal 
forces  shown  in  Figures  4.2  and  4.6.  As  can  be  seen  in  Figures  4.2a  and  4.6a,  the 
vertical  excavation  forces  from  the  method  of  force  residuals  and  the  Clough  and 
Duncan  procedure  are  the  same  for  the  first  stage  of  excavation.  Also,  the  sum¬ 
mation  of  the  vertical  excavation  forces  equals  24,000  lb  for  both  procedures, 
which  is  the  same  as  the  weight  of  the  two  excavated  elements.  Thus,  force 
equilibrium  is  satisfied  for  both  procedures  for  the  first-stage  excavation. 
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Figure  4.8.  Interpolated  stresses  at  element  edges  (a)  and  equivalent  nodal  forces  (b)  along  the  excava¬ 
tion  boundaries  from  Element  2,  calculated  from  the  Clough  and  Duncan  procedure  (stresses 
are  in  kilopascals,  forces  are  in  kilonewtons) 


=  lOkPa 


1  m 


(a)  (b) 


Figure  4.9.  Nodal  forces  required  to  satisfy  equilibrium  offerees  in  the  (a)  vertical  and  (b)  horizontal 
directions 
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The  lack  of  force  equi¬ 
librium  for  the  Clough  and 
Duncan  procedure  becomes 
apparent  for  the  second  stage 
of  excavation.  This  can  be 
seen  by  comparing  the  verti¬ 
cal  excavation  forces  in 
Figures  4.2b  and  4.6b.  The 
summation  of  the  vertical 
excavation  forces  gives 
24,000  lb  (=  6,324  +  14,347 
+  7,647 -4,3 18)  for  the 
method  of  force  residuals. 

This  value  is  again  equal  to  Figure  4. 1 0.  Equivalent  nodal  forces  based  on 

the  total  weight  of  the  two  minimization  of  force  residuals 

excavated  elements  in  the  (forces  are  in  kilonewtons) 

second  stage  of  excavation. 

In  comparison,  the  summation  of  the  vertical  forces  from  the  Clough  and  Duncan 
procedure  is  equal  to  30,397  lb  (=  6,800  +  15,030  +  9,555  -  988).  This  is  much 
larger  than  the  total  weight  of  the  two  excavated  elements  of  24,000  lb  and, 
therefore,  force  equilibrium  is  not  satisfied  for  the  second  stage  of  excavation. 


4.2.3  Convergence  of  solution 

As  pointed  out  by  Ghaboussi  and  Pecknold  (1984),  a  clear  distinction  should 
be  made  between  solution  accuracy  and  step-size  independence.  While  it  is  true 
that  more  accurate  results  should  be  obtained  by  the  use  of  higher  order  elements 
or  finer  discretizations,  this  has  little  to  do  with  the  problem  of  cumulative  step 
errors.  If  the  numerical  solution  process  is  carried  out  correctly,  the  solution  will 
always  be  independent  of  the  number  of  steps,  no  matter  how  crude  the  finite 
element  mesh  used  to  model  the  stress  distribution.  Of  course,  a  finer  discretiza¬ 
tion  will  yield  a  more  accurate  result  using  the  method  of  force  residuals.  How¬ 
ever,  as  pointed  out  by  Brown  and  Booker  (1986)  the  method  of  force  residuals 
does  not  require  a  large  degree  of  mesh  refinement  or  the  use  of  high-order  finite 
elements  to  obtain  accurate  results  for  an  elastic  material.  This  was  shown  above 
where  the  finite  element  model  using  four-noded  elements  yielded  results  almost 
identical  to  the  finite  element  model  using  eight-noded  finite  elements  (and 
doubling  the  number  of  nodes). 


4.3  Implementation  Issues  in  SOILSTRUCT- 
ALPHA 

Based  on  the  comparisons  made  in  Section  4.1  and  the  results  of  the  finite 
element  analysis  carried  out  in  Section  4.2,  it  is  recommended  that  the  method  of 
force  residuals  be  implemented  as  another  option  for  excavation  modeling  in 
SOILSTRUCT-ALPHA.  There  are  several  alternatives  for  implementing  the 
method  of  force  residuals  in  SOILSTRUCT-ALPHA.  In  the  different  alterna¬ 
tives,  the  overriding  issue  is  the  type  of  element  that  is  available  or  can  be 
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implemented  in  SOILSTRUCT-ALPHA.  Currently,  the  only  continuum  element 
available  in  this  code  is  the  so-called  QM5  quadrilateral  element.  The  properties 
of  this  element  are  discussed  further  below. 

The  alternatives  for  implementing  the  method  of  force  residuals  in 
SOILSTRUCT-ALPHA  are 

a.  Using  the  existing  QM5  finite  element  routine  with  constant  element 
stresses. 

b.  Using  the  existing  QM5  finite  element  routine  but  with  element  stresses 
calculated  at  the  Gaussian  quadrature  points. 

c.  Using  a  four-noded  isoparametric  finite  element. 

d.  Using  higher  order  (e.g.,  eight-noded)  finite  elements. 

The  simplest  option  is,  of  course,  to  retain  the  existing  Clough  and  Duncan 
procedure  in  SOILSTRUCT-ALPHA,  and  use  as  few  excavation  steps  as  pos¬ 
sible  based  on  the  recommendation  by  Christian  and  Wong  (1973).  As  shown 
above,  a  single-stage  excavation  may  yield  results  comparable  to  those  obtained 
from  the  method  of  force  residuals,  at  least  for  linear  elastic  problems.  Because 
this  approach  is  not  feasible  for  many  problems  of  practical  interest,  it  is  recom¬ 
mended  that  the  method  of  force  residuals  be  implemented  in  SOILSTRUCT- 
ALPHA.  The  following  sections  discuss  issues  related  to  this  implementation. 


4.3.1  Implementation  using  QM5  or  other  element  types 

Using  the  existing  QM5  element  appears  to  be  the  most  straightforward  way 
of  implementing  the  method  of  residuals  in  SOILSTRUCT-ALPHA.  It  is  desir¬ 
able  to  retain  the  use  of  QM5  element  since  this  would  involve  the  least  amount 
of  modifications  to  SOILSTRUCT-ALPHA,  and  there  is  already  a  large  experi¬ 
ence  base  with  the  use  of  this  element.  In  the  following,  a  brief  description  of 
QM5  is  given  first  to  clarify  the  feasibility  of  implementing  the  method  of  force 
residuals  in  SOILSTRUCT-ALPHA. 


QM5  is  one  of  the  family  of  higher  order  quadrilateral  elements  developed 
by  Doherty,  Wilson,  and  Taylor  (1969)  for  stress  analysis  of  axisymmetric  solids. 
The  axisymmetric  formulation  was  reduced  to  a  plane  strain  formulation  when 
QM5  was  implemented  in  SOILSTRUCT.  The  main  feature  of  the  five-noded 
quadrilateral  element  QM5  is  that  it  has  better  bending  characteristics  over  other 
elements  with  the  same  number  of  nodes.  This  is  accomplished  by  using  a  dis¬ 
placement  function  of  the  following  form: 


u  =  a0+als  +  a2t  +  a3st  +  a4(l-s2)(l-t 2) 
v  =  ba  +  bts  +  b2t  +  b3st  +  b4  (l  -  s2 )  ( 1  - 12 ) 


(4.1) 
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where  u  and  v  are  the  horizontal  and  vertical  displacements,  a0  ...  at,  and  b„  ...  b4 
are  interpolation  constants  that  can  be  determined  from  the  values  of  u  and  v  at 
the  five  nodes  (the  fifth  node  is  located  at  the  element  centroid),  and  s  and  t  are 
the  local  element  coordinates.  (See  Appendix  A  for  an  explanation  of  local 
coordinates.)  The  first  four  terms  in  the  two  polynomials  in  Equation  4. 1  are  the 
same  as  those  used  in  standard  four-noded  isoparametric  elements.  Note  that  the 
equation  is  an  incomplete  polynomial. 

The  fifth  terms  involving  a4  and  b4  are  added  to  improve  the  bending  per¬ 
formance  of  the  standard  four-noded  quadrilateral  element.  However,  compati¬ 
bility  requires  that  only  the  first  four  bilinear  terms  in  Equation  4.1  can  be  used 
for  the  four  comer  nodes.  This  forces  the  element  to  have  linear  displacements 
along  the  element  sides.  To  retain  the  improved  bending  performance  of  the 
element,  a  constant  shear  strain  is  imposed  on  the  element.  The  shear  strain  is 
calculated  at  the  center  or  the  fifth  node  of  the  element  using  the  displacement 
function  given  in  Equation  4.1.  Once  the  shear  strain  has  been  calculated,  the 
central  node  is  deleted  from  the  finite  element  formulation  using  a  method  of 
static  condensation.  In  the  assembly  of  the  stiffness  matrix,  only  the  four  comer 
nodes  are  included,  and  the  central  node  does  not  participate  in  the  displacement 
calculations. 

Due  to  the  assumed  constant  shear  strain  in  the  element,  only  one  set  of 
stresses  is  currently  calculated  in  QM5,  and  these  are  calculated  at  the  element 
center.  This  is  equivalent  to  assuming  that  QM5  is  a  constant  stress  element.  It  is 
possible  to  calculate  at  least  four  different  stresses  at  four  Gaussian  stress  points 
in  QM5;  however,  this  will  require  the  displacements  at  the  central  node  to  be 
calculated.  Two  methods  were  tried  to  recover  the  displacements  at  the  central 
node  from  the  four  comer  nodes:  averaging  the  four  comer  displacements,  and 
using  a  technique  to  decondensate  the  internal  degrees  of  freedom  (Desai  and 
Abel  1972). 

The  first  method,  while  very  simple  and  straightforward,  yielded  very  high 
forces  at  the  central  node  that  are  an  order  of  magnitude  larger  than  the  forces  at 
the  comers. 

The  second  method  was  found  to  be  very  cumbersome,  involving  re-solution 
of  the  systems  of  equations  in  order  to  recover  the  internal  degrees-of-freedom 
from  the  four  comer  nodes.  Such  a  procedure  was  nonetheless  tried,  and  the 
resulting  finite  element  formulation  was  used  in  the  same  problem  analyzed 
(shown  in  Figure  4.2).  To  achieve  complete  integration  of  the  fourth-order  dis¬ 
placement  function  (Equation  4.1),  a  3  by  3  Gaussian  quadrature  was  used  to 
integrate  the  stiffness  matrix,  and  the  element  stresses  are  determined  at  the  nine 
quadrature  points.  The  result  of  a  two-stage  excavation  using  the  QM5  with 
recovered  internal  displacements  is  shown  in  Figure  4.1 1.  As  can  be  seen,  the 
overall  pattern  of  displacement  looks  reasonable,  but  individual  elements  gave 
very  unusual  displacement  modes.  The  displacement  mode  of  each  element  is 
due  to  the  imposed  constant  shear  strain  within  the  element  and  the  compatibility 
condition  that  elements  must  show  linear  displacements  along  the  edges.  The 
results  of  the  analysis  indicate  that  it  is  difficult  to  obtain  a  reasonable  displace¬ 
ment  response  from  the  QM5  element  if  the  internal  displacements  are  recovered 
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Figure  4.11.  Deformed  mesh  after  excavation  using  the  QM5  element  with 
recovered  internal  displacements 

to  calculate  stresses  at  the  Gaussian  points  while  at  the  same  time  keeping  the 
constant  shear  strain  condition. 

Regardless  of  how  the  internal  element  stresses  are  determined,  only  one  set 
of  element  nodal  forces  can  be  calculated  for  an  element.  This  means  that, 
although  stresses  can  be  calculated  at  four  or  nine  Gaussian  points,  these  stresses 
will  have  to  be  weighted  and  summed  to  get  a  single  set  of  nodal  forces.  How¬ 
ever,  calculating  stresses  at  the  Gaussian  points  will  give  a  better  representation 
of  the  stresses  within  the  elements,  and  thus  highly  stressed  regions  in  the  ele¬ 
ment  are  better  represented  than  by  assuming  a  constant  stress  distribution.  A 
better  representation  of  element  stresses  will  eventually  yield  a  more  realistic 
distribution  of  element  forces.  Finally,  it  should  be  noted  that  the  central  node 
should  give  zero  residual  forces  at  the  end  of  the  gravity  tum-on  and  excavation 
analyses  since  a  central  node  is  not  connected  to  any  other  node.  Thus,  there  is  no 
added  value  in  keeping  the  internal  degree-of-freedom  in  the  calculation  of 
excavation  forces. 

Based  on  the  above  observations  and  considerations  of  other  element  types,  a 
summary  of  the  finite  element  options  for  simulating  excavation  in 
SOILSTRUCT- ALPHA  is  given  in  Table  4.3 . 
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Table  4.3 

Summary  of  Advantages  and  Disadvantages  of  Finite  Element 
Routines  That  Could  Be  Used  in  Implementing  the  Method  of  Force 
Residuals  in  SOILSTRUCT-ALPHA 

Method 

Advantages/Disadvantages 

Using  QMS  with  only 
one  set  of  element 
stresses. 

This  is  the  simplest  and  most  straightforward  approach  of  introducing  the 
method  of  force  residuals  in  SOILSTRUCT-ALPHA. 

Most  of  the  required  routines  are  already  available  in  the  code.  j 

It  is  anticipated  that  even  with  a  constant-stress  QM5  element,  there  will 
be  improvement  over  the  Clough  and  Duncan  procedure. 

Using  QM5  but  extract 
four  stresses  at  the  four 
Gaussian  points. 

This  will  require  very  complicated  programming  to  recover  the  internal 
degree-of-freedom  (DOF)  that  was  taken  out  during  the  static 
condensation. 

The  extra  calculations  will  significantly  slow  down  the  code. 

The  added  sophistication  is  shown  not  to  give  substantial  improvement 
from  the  current  constant  stress  QMS.  In  fact,  the  QMS  element  with 
recovered  internal  displacements  is  shown  to  exhibit  rather  unusual 
displacement  response. 

Introduce  a  new  four- 
noded  element 

This  is  preferred  over  the  second  option,  as  it  will  require  less 
programming,  and  the  code  will  be  more  robust  than  one  where  internal 
DOFs  are  recovered  from  QM5.  This  is  particularly  true  if  elements  with 
four  nodes  are  used  as  in  the  reduced  QM5  element. 

An  isoparametric  four-noded  quadrilateral  (Q4),  which  will  yield  four 
stresses  at  the  Gaussian  points,  can  replace  or  be  used  side-by-side  with 
the  QMS  that  has  the  same  number  of  DOFs. 

The  Q4  does  not  have  the  advantage  of  the  improved  bending 
performance  of  the  QMS.  The  Q4  can  be  modified  to  produce  an 
improved  bending  performance,  again  by  imposing  a  constant  strain 
condition  in  the  element,  as  was  done  in  the  element  QM4  developed  by 
Doherty,  Wilson,  and  Taylor  (1969).  As  in  the  QMS,  the  penalty  is  that  only 
one  set  of  stresses  can  be  calculated  in  the  QM4. 

Introduce  a  higher  order 
element. 

Higher  order  elements  will  yield  better  performance  including  bending  that 
lower  order  elements. 

This  will  require  the  most  modification,  particularly  in  the  data  structure  of 
the  code  to  allow  for  midside  nodes.  Modifications  would  be  necessary  in 
almost  all  parts  of  the  code,  including  the  interface  element  routine  and 
the  calculation  of  initial  stresses 

4.3.2  Interface  elements 

Another  issue  that  needs  to  be  considered  in  implementing  a  new  excavation 
routine  in  SOILSTRUCT-ALPHA  is  the  influence  of  interface  elements.  These 
elements  are  introduced  between  the  soil  and  a  retaining  structure  to  permit  rela¬ 
tive  deformation  between  the  soil  and  the  structure.  There  has  been  an  active 
development  of  models  for  interface  behavior  in  SOILSTRUCT-ALPHA 
(Gomez,  Filz,  and  Ebeling  2000a,  2000b).  Currently,  however,  the  interface 
element  has  no  influence  on  the  excavation  routine  in  SOILSTRUCT-ALPHA. 
During  excavation,  the  interface  elements  are  assumed  to  be  “locked”  and  are 
given  large  normal  and  shear  stiffness  values  during  the  excavation  simulation. 

Interface  elements  are  1-D  elements  or  elements  without  thickness  intro¬ 
duced  between  the  structure  and  the  soil.  Interface  elements  allow  adjacent 
continuum  soil  and  structure  elements  to  move  independently  of  each  other.  In 
the  case  of  the  problem  shown  in  Figures  1.3  and  2.1,  interface  elements  would 
be  used  between  the  two  sides  of  the  slurry  wall  and  the  soil  (see  Figure  4.12). 
The  behavior  of  these  interface  elements  is  governed  by  normal  and  shear 
stiffnesses,  with  the  former  controlling  the  opening  and  closure  and  the  latter 
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Figure  4.12.  Use  of  interface  elements  to  model  interaction  between  soil  and 
structure 

controlling  the  shear  displacement  along  the  interface.  In  SOILSTRUCT- 
ALPHA,  the  interface  element  has  four  nodes,  with  each  of  two  nodes  in  a  pair 
having  the  same  coordinates. 

Interface  behavior  can  have  important  influence  on  the  excavation  forces  that 
are  built  up  in  the  structure  due  to  the  backfilling  of  the  soil  around  the  structure, 
or  due  to  the  driving  shear  stresses  in  the  case  of  structures  driven  into  the  soils. 
In  the  case  of  backfilling,  an  important  effect  is  the  generation  of  downdrag  or 
vertical  shear  force  exerted  by  the  backfill  on  the  structure,  as  referred  to  in 
Section  1 .2.  As  noted  above,  downdrag  effects  have  been  reported  in  earlier 
studies  by  Ebeling,  Duncan,  and  Clough  (1990);  Ebeling  et  al.  (1993);  Ebeling 
and  Mosher  (1996);  Ebeling  and  Wahl  (1977);  and  Ebeling,  Peters,  and  Mosher 
(1997).  The  main  purpose  of  interface  elements  is  to  be  able  to  permit  slippage 
between  the  soil  and  the  structure  during  construction  and  also  during  subsequent 
loading  of  the  structure. 

To  include  interface  elements  in  the  excavation  simulation  using  the  method 
of  force  residuals,  it  is  necessary  to  calculate  the  internal  forces  in  the  interface 
elements.  In  case  of  continuum  elements,  the  internal  forces  are  calculated  from 
the  built-up  stresses  {a}  in  the  element  using  the  equation 

Rt}=  dV  (4.2) 


A  more  direct  approach  can  also  be  used  to  calculate  the  internal  force  that  is 
applicable  to  both  continuum  and  interface  elements.  Again,  for  continuum 
elements,  the  stress  increments  {a}  are  related  to  the  strain  increments  {e}  by  the 
constitutive  matrix  [D\. 
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{cHl>]{e} 


(4.3) 


while  the  strain  increments  {£}  are  related  to  the  incremental  element  nodal 
displacements  {*/}  by  the  strain-displacement  matrix: 

{e}  =  [B]{»'}  (4.4) 

Substitution  in  Equation  4.2  gives  the  internal  nodal  loads: 

[F„}=  \{B}T  [D]{B}{u-}dV  (4.5) 

V 

Here,  the  quantity  J{i?}r  [D\{B\dV  is  actually  equal  to  the  element  stiffness 

v 

matrix  \jce  ]  and,  therefore.  Equation  4.4  can  also  be  written  as 


W-MM’  [*']=/ W[D\{B}dV  (4.6) 

V 

Equation  4.6  provides  another  way  of  calculating  equivalent  internal  nodal 
forces  from  element  stresses.  In  the  case  of  interface  elements,  the  internal  forces 
can  be  calculated  from  the  relative  nodal  displacements  along  the  element  and  the 
normal  and  shear  stiffnesses  of  the  interface: 

<4-7> 

Note  that  Equation  4.7  is  applicable  to  all  kinds  of  elements,  and  thus,  the 
internal  forces  can  be  computed  in  a  consistent  manner  and  included  in  the  inter¬ 
nal  stresses  from  the  continuum  elements.  The  internal  forces  in  the  interface  are 
then  included  in  the  calculation  of  the  residual  forces  due  to  excavation.  Deletion 
of  interface  elements  causes  no  change  in  the  external  body  loads  as  interface 
elements  have  no  self-weights.  In  this  manner,  no  special  programming  is  needed 
to  handle  interface  elements  in  the  excavation  simulation  using  the  method  of 
force  residuals.  In  the  case  of  nonlinear  material  response,  Equation  4.7  must  be 

evaluated  incrementally  using  the  nonlinear  stiffness  matrix  J  consistent  with 
the  variation  of  stresses  and  strain. 

To  illustrate  the  importance  of  interface  response  and  to  show  the  use  of  the 
method  of  force  residuals  to  account  for  interface  elements  in  excavations,  the 
problem  shown  in  Figure  4.13  below  is  used.  This  is  a  similar  problem  to  that 
analyzed  in  Section  4.3.1  except  that  interface  elements  are  placed  along  the  ver¬ 
tical  excavation  boundaiy,  extending  from  the  top  to  the  bottom  of  the  model.  In 
practice,  such  elements  will  be  installed  between  soil  elements  and  a  structural 
element  (e.g.,  a  slurry  wall).  However,  in  this  simple  problem,  the  interface  ele¬ 
ments  are  placed  artificially  between  soil  elements.  Again,  the  objective  is  to 
simulate  the  deformation  of  the  soil  after  Elements  1,  5, 2,  and  6  are  deleted. 
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The  interface  element  used  is  the  four-noded  isoparametric  interface  element 
developed  by  Beer  (1985)  and  Carol  and  Alonso  (1983).  This  interface  element  is 
similar  in  formulation  to  the  interface  element  developed  by  Goodman,  Taylor, 
and  Brekke  (1968)  with  the  numerical  integration  proposed  by  Morrison  and 
Duncan  (1995).  For  simplicity,  the  relative  displacements  at  the  two  sides  of  the 
element  are  related  to  the  normal  and  shear  forces  Fn  and  Fs  along  the  element 
by  linear  elastic  normal  and  shear  stiffnesses  Kn  and  Ks : 


£ 

1 

X  o' 

XJ 

.0  ks 

where  [ulop  -  ubM )  and  ( vlnp  -  vbnl )  are  the  relative  displacements  normal  and 
parallel  to  the  interface. 

The  response  of  the  soil  after  the  single-stage  excavation  is  shown  in 
Figures  4.14-4.16.  In  all  the  simulations,  the  same  normal  stiffness  of 
Kn  =  1-106  lb/ft  is  used,  while  three  different  values  of  the  shear  stiffness 

Ks  (equal  to  MO6,  MO3,  11 02  lb/ft)  were  used.  The  values  of  the  maximum 

vertical  displacements  or  heave  at  the  base  of  the  excavation  from  the  simulations 
with  the  different  shear  stiffness  values  are  also  summarized  in  Table  4.4  below. 
The  simulation  using  Ks  =  1  •  1 06  lb/ft  gave  almost  the  same  response  as  the  case 

without  interface  elements,  except  for  the  slight  interpenetration  or  overlap  of  the 
soil  elements  attached  to  the  interface  below  the  excavation.  Note  that  the 
overlap  of  the  finite  elements  shown  in  Figure  4.14  can  be  reduced  by  using  a 
higher  value  of  the  contact  normal  stiffness  Kn . 


Table  4.4 

Effects  of  Interface  Shear  Stiffness  on  Maximum  Heave  at  Base  of 
Excavation 

1  Shear  Stiffness 

Maximum  Vertical  Displacement 

J  Ks  =1-10°  lb/ft 

0.403  ft 

A>1103  lb/ft 

0.402  ft 

AT,  =1-102  lb/ft 

0.398  ft 

|  Notes:  Kn  =  110B 

lb/ft  for  al!  simulations.  To  convert  feet  to  meters,  multiply  by  0.3048.  | 

As  shown  in  Table  4.4,  there  is  very  little  effect  of  the  interface  shear  stiff¬ 
ness  on  the  maximum  heave  at  the  base  of  the  excavation.  The  main  effect  of  the 
shear  stiffness  is  on  the  pattern  of  deformation  below  the  excavation.  The  heave 
of  the  remaining  soil  below  the  excavation  becomes  more  uniform  as  the  inter¬ 
face  shear  stiffness  is  reduced.  For  the  lowest  value  of  Ks  =  11 02  lb/ft  used,  the 

soil  below  the  excavation  deformed  in  a  piston-like  manner  with  an  almost  com¬ 
plete  vertical  slip  along  the  interface. 
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The  results  shown  in  Figures  4.13-4.16  show  the  significance  of  including 
interface  elements  in  the  excavation  simulation.  The  results  also  demonstrate  that 
interface  elements  can  be  included  in  a  consistent  manner  using  the  method  of 
force  residuals  with  internal  forces  determined  from  Equation  4.7. 


4.3.3  Nonlinear  material  properties  and  the  ALPHA  method 

Soil  and  interface  behavior  are  modeled  using  nonlinear  models  in 
SOILSTRUCT-ALPHA.  For  instance,  soils  are  modeled  using  the  Duncan  and 
Chang  (1970)  hyperbolic  model.  Currently,  the  so-called  ALPHA  method 
(Ebeling,  Duncan,  and  Clough  1990;  Regalado,  Duncan,  and  Clough  1992)  is 
used  in  the  code  for  analyzing  the  nonlinear  response  of  soil  structures.  It  is, 
therefore,  important  to  investigate  how  the  method  of  force  residuals  can  be 
implemented  in  SOILSTRUCT-ALPHA  taking  into  consideration  the  ALPHA 
method.  Nonlinear  analysis  is  carried  out  in  a  series  of  increments  or  steps  in 
SOILSTRUCT-ALPHA.  For  each  increment  the  deformation  moduli  of  the 
different  materials  are  adjusted  according  to  the  level  of  the  current  stress  and 
deformation,  and  the  analysis  for  each  increment  is  carried  out  as  if  the  materials 
behave  linearly. 

Consider  the  case  of  the  structure  shown  in  Figure  4.17.  Let  P  be  the  current 
load,  let  A P  be  the  load  increment,  and  assume  that  the  current  stresses  at  the 
most  stressed  element  are  below  the  failure  surface  in  the  Mohr-Coulomb 
diagram  (Figure  4.17a).  Due  to  the  finite  nature  of  the  load  increment  A P,  some 
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Figure  4. 1 6.  Deformed  mesh  after  excavation  with  K,  =  1  •  1 02  Ib/ft  (deformations 
magnified  lOx) 

elements  will  experience  stresses  that  overshoot  the  failure  surface  (Fig¬ 
ure  4.17b).  Ideally,  the  load  should  be  applied  in  each  increment  such  that  it  is 
large  enough  to  bring  the  most  severely  stressed  soil  element(s)  exactly  to  failure 
(i.e.,  no  overshoot),  as  shown  in  Figure  4.17c.  The  load  A P  should  then  be 
reduced  by  a  factor  in  order  to  bring  the  most-stressed  element  just  to  a  failure 
condition.  This  factor,  called  a  (hence  the  term  “ALPHA  method”),  is  calculated 
from  the  stresses  that  are  just  tangential  to  the  Mohr-Coulomb  surface  (Fig¬ 
ure  4.17c).  Using  the  a-factor  to  reduce  the  load  will  ensure  that  failure  takes 
place  only  in  the  worst  stressed  element  and  nowhere  else  in  the  modeled  region. 
The  ALPHA-factor  is  applied  to  both  continuum  and  interface  elements. 

Due  to  its  incremental  nature,  the  method  of  force  residuals  can  be  equally 
and  directly  implemented  using  the  ALPHA-method  to  integrate  the  nonlinear 
equations.  As  shown  in  Section  3.2,  in  the  method  of  force  residuals,  the  loads 
due  to  excavation  are  calculated  from  the  difference  in  the  internal  and  external 
forces: 


[*]{*,}={*} 

(4.9) 

{R}={Fmt}-{Fexl} 

(4.10) 
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Figure  4.17.  ALPHA-method  in  factoring  out  load  increments  to  avoid  overshoot 


Therefore,  using  the  ALPHA-method,  the  a-factor  can  be  directly  applied  to  the 
residual  force  vector  {/?} .  An  advantage  is  that  the  force  residuals  may  either  be 

applied  using  a  single  stage  excavation  or  using  multiple  stages  of  excavation.  In 
a  single-stage  excavation,  the  residual  forces  along  the  excavation  boundary  are 
calculated  all  at  once  but  are  applied  incrementally  using  the  a-factor.  An  alter¬ 
native  is  to  perform  the  excavation  in  stages  taking  a  few  elements  at  a  time  (via 
a  layer-by-layer  excavation  for  instance),  and  apply  the  ALPHA-method  for  each 
excavation  stage.  The  latter  method  is  preferred,  as  the  excavation  response  of 
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the  soil  will  now  be  dependent  on  the  number  of  excavation  stages  used  and  the 
size  of  each  excavation  stage. 


4.3.4  Programming  issues 

There  are  some  programming  issues  that  need  to  be  addressed  in  the  imple¬ 
mentation  of  the  method  of  force  residuals  in  SOILSTRUCT-ALPHA.  As 
mentioned  above,  implementing  the  method  of  residuals  using  the  existing  QM5 
element  with  constant  stresses  would  entail  the  least  amount  of  modification  of 
SOILSTRUCT-ALPHA.  The  modifications  and  related  programming  issues 
discussed  below  are  based  on  using  the  QM5  with  constant  stresses. 

The  Clough  and  Duncan  excavation  procedure  is  currently  implemented  in 
the  subroutine  EXCAV.  This  subroutine  is  called  once  every  loading  increment 
by  the  MAIN  program.  For  each  element  to  be  excavated,  EXCAV  locates  the 
centers  of  the  four  interpolation  elements  and  performs  the  interpolation  on  the 
stresses  at  the  centers  of  these  elements.  The  interpolation  constants  are  then  used 
to  determine  the  stresses  at  the  nodes  of  the  elements,  which  are  then  passed  to 
the  subroutine  EQNDFO,  which  converts  the  nodal  stresses  to  equivalent  nodal 
forces  along  the  excavation  boundary. 

Instead  of  rewriting  EXCAV,  it  is  proposed  to  develop  a  new  subroutine 
EXCAV2  to  implement  the  method  force  residuals  in  SOILSTRUCT-ALPA. 
EXCAV  can  be  retained,  and  EXCAV2  may  be  used  as  an  alternative  to  EXCAV 
so  that  users  will  be  allowed  to  choose  either  excavation  procedure.  Another 
solution  is  to  develop  a  new  version  of  SOILSTRUCT-ALPHA  where  EXCAV2 
completely  replaces  EXCAV  and  is  the  only  procedure  available. 

The  main  calculations  to  be  done  in  EXCAV2  are  these:  (a)  modifying  the 
element  properties  of  the  excavated  elements,  (b)  calculating  the  internal  forces, 
(c)  modifying  the  load  vector  to  account  for  deleted  elements,  and  (d)  calculating 
the  residual  force  vector.  The  calculation  of  the  internal  element  forces  will  be 
performed  after  the  incremental  displacements  have  been  solved  for  a  previous 
loading  increment  and  the  global  displacements  and  element  stresses  have  been 
updated.  The  calculations  of  the  residual  forces  should  be  made  only  on  the 
nodes  along  the  excavation  boundary.  The  residual  forces  are  incremental,  and 
the  calculated  displacements  and  stresses  are  also  incremental  and  must  be  added 
to  the  cumulative  displacement  and  stress  vectors.  Special  precautions  must  be 
made  to  incorporate  the  effects  of  seepage,  temperature,  and  added  elements 
from  fill  placement  lifts  in  the  residual  force  calculations.  The  flow  charts  for  two 
recommended  procedures  to  incorporate  the  method  of  residuals  in 
SOILSTRUCT-ALPHA  are  given  in  Figures  4.18  and  4.19. 

4.3.4.1  Calculation  of  equivalent  internal  nodal  forces.  There  are  two 
options  in  calculating  the  internal  element  force  vector:  (a)  using  the  [J5]  -matrix, 

the  element  stresses  {a} ,  and  Equation  3.2  (as  illustrated  in  Figure  4.18)  and 
(b)  using  the  element  matrices  [ke],  the  element  displacements {we} ,  and 
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Figure  4. 1 8.  Flowchart  for  implementing  the  method  of  force  residuals  in  SOILSTRUCT- 
ALPHA  using  element  stresses  for  calculating  internal  element  forces 


78 
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Figure  4.19.  Flowchart  for  implementing  the  method  offeree  residuals  in  SOILSTRUCT- 

ALPHA  using  incremental  element  displacements  in  calculating  internal  element 
forces 
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Equation  4.7  (as  illustrated  in  Figure  4.19).  The  element  stiffness  matrices  are 
calculated  in  the  subroutines  QUAD,  JTSTF,  and  BAREL  for  the  continuum, 
interface,  and  structural  elements,  respectively.  These  subroutines  pass  the 
stiffness  matrices  to  the  subroutine  STRSTF,  which  assembles  the  global  stiff¬ 
ness  matrix  at  the  start  of  every  new  load  increment.  Although  QUAD,  JTSTF, 
and  BAREL  generate  the  [Are]  -matrices,  the  element  [5]  -matrices  are  not  gen¬ 
erated  explicitly  (QUAD  assembles  the  elements  of  the  [5]  -matrix,  but  these  are 

not  stored  in  any  matrix  and  are  currently  not  passed  to  other  subroutines).  Ele¬ 
ment  stresses  are  generated  in  STRESS,  JSTRES,  and  BSTRES  for  the  con¬ 
tinuum,  interface,  and  structural  elements,  respectively.  Note  that  [jte]  is  the 
tangent  element  stiffness  matrix  for  incremental  analyses. 

Using  the  first  option  would,  therefore,  require  that  the  [fi] -matrices  be 

generated  (for  JTSTF  and  BAREL),  stored  in  a  matrix,  and  passed  via  the 
COMMON  block.  In  the  second  option,  the  internal  forces  must  be  stored  in  an 
array  and  updated  incrementally  using  the  incremental  displacements  from  each 
increment  of  loading  (and  taking  into  consideration  the  nonlinear  response  of  the 
soil  and  the  interface  elements).  The  accumulated  internal  forces  should  include 
the  internal  forces  from  displacements  generated  in  the  gravity  turn-on  analysis. 

The  advantage  of  calculating  the  internal  element  forces  by  using  the  element 
stresses  and  Equation  3.2  is  that  the  element  stresses  are  already  updated  and 
stored  globally  by  SOILSTRUCT- ALPHA.  The  disadvantage  of  the  approach  is 
that  the  element  [5]  matrices  must  be  generated;  however,  these  are  needed  only 
for  the  excavated  elements. 

The  advantage  of  calculating  the  internal  element  forces  by  accumulating  and 
storing  the  element  internal  forces  from  Equation  4.7  using  the  element  tangent 
stiffness  matrix  [ke  ]  is  that  a  new  matrix  corresponding  to  [fi]  is  not  created  and 

passed  via  the  COMMON  statement.  The  disadvantage  of  this  approach  is  that 
the  internal  element  load  vector  must  be  stored  for  each  element  in  the  entire 
mesh  and  updated  continuously. 

43.4.2  Calculation  of  external  element  forces.  To  calculate  the  external 
forces,  a  new  subroutine  must  be  written  that  modifies  the  current  external  load 
vector  to  account  for  the  deletion  of  the  weight  of  the  excavated  elements. 
Changes  in  the  external  loads  due  to  deletion  of  the  weight  of  structural  elements 
can  and  should  also  be  accounted  for.  This  requires  that  the  load  vector  must  be 
stored  separately  and  not  be  overwritten  by  the  displacement  vector  as  is  com¬ 
monly  done  in  many  finite  element  programs.  The  load  vector  must  also  continu¬ 
ously  track  the  changes  made  to  the  external  loads  calculated  from  the  in  situ 
stresses  created  in  the  INITIAL  subroutine.  These  changes,  which  are  mainly 
from  the  weights  of  added  elements  simulating  embankment  lift,  are  calculated  in 
the  subroutine  BUILD.  Changes  from  external  loads  on  the  boundaries  of  the 
excavated  elements  should  also  be  accounted  for  in  the  load  vector. 
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5  Summary  and  Conclusions 


The  Corps  uses  SOILSTRUCT-ALPHA  to  perform  soil-structure  interaction 
analyses  of  multi-anchored  or  tieback  retaining  walls.  Permanent  multi-anchored 
walls  have  been  used  as  guide  walls  and  approach  walls  on  navigation  projects, 
and  as  retaining  walls  on  highway  and  railroad  protection  and  relocation  projects. 
A  critical  aspect  for  SSI  analyses  of  such  structures  is  the  modeling  of  the 
excavation  process.  The  case  of  the  Bonneville  temporary  tieback  wall  was  used 
to  illustrate  the  importance  of  SSI  analyses,  and  particularly  the  need  for  a  robust 
and  efficient  procedure  to  simulate  excavation  processes.  The  excavation 
algorithm  in  SOILSTRUCT  was  developed  in  the  late  1960s  by  Clough  and 
Duncan  (1969),  and  it  has  not  been  updated  since  that  time,  even  though 
improved  methods  have  been  developed  by  Ghaboussi  and  Pecknold  (1984), 
Borja  et  al.  (1989),  and  others. 

The  objectives  of  the  study  described  in  this  report  are  to 

a.  Complete  a  comprehensive  literature  review  of  numerical  modeling  of 
excavation. 

b.  Determine  the  suitability  of  existing  excavation  algorithms  for  their  use 
in  SOILSTRUCT. 

c.  Provide  the  information  necessary  to  implement  an  improved  excavation 
algorithm  in  SOILSTRUCT. 

The  Clough  and  Duncan  excavation  algorithm  was  reviewed  extensively.  An 
example  problem  was  analyzed  to  illustrate  the  main  steps  in  this  algorithm, 
particularly  the  extrapolation  of  stresses  from  nearby  elements  to  the  excavation 
boundary  and  integration  of  the  extrapolated  stresses  to  produce  excavation 
forces  at  element  nodes.  In  addition  to  the  Clough  and  Duncan  procedure,  other 
algorithms  were  also  reviewed,  including  1)  the  method  of  force  residuals,  2)  the 
method  based  on  accumulation  of  excavation  forces  along  the  excavation  bound¬ 
aries,  and  3)  the  method  based  on  hybrid  finite  elements.  Of  these  procedures,  the 
method  of  force  residuals  is  the  most  widely  used  procedure  for  modeling  exca¬ 
vation.  The  example  problem  analyzed  using  the  Clough  and  Duncan  procedure 
was  reanalyzed  to  illustrate  the  main  steps  involved  in  the  method  of  force 
residuals. 

The  different  excavation  algorithms  were  rigorously  analyzed  and  compared 
in  terms  of  their  numerical  performance.  It  was  shown  that  the  Clough  and 
Duncan  excavation  procedure  yields  results  that  are  dependent  on  the  number  of 
excavation  stages  used  in  the  simulations  for  linear  elastic  materials.  For  such 


Chapter  5  Summary  and  Conclusions 


materials,  it  has  been  shown  formally  by  Ishihara  (1970)  that  the  results  of  the 
excavation  should  be  independent  of  the  number  of  steps  and  the  sequence  used 
to  simulate  the  excavation.  The  lack  of  uniqueness  in  the  results  of  the  Clough 
and  Duncan  excavation  algorithms  stems  from  the  fact  that  the  excavation  forces 
applied  along  the  excavation  boundary  are  not  always  consistent  with  the  stresses 
in  the  excavated  elements.  Force  equilibrium  is  not  necessarily  satisfied,  and  the 
resulting  deformations  can  be  either  larger  or  smaller  than  they  should  be.  The 
lack  of  force  equilibrium  can  result  in  a  distribution  of  excavation  unloading 
force  in  which  either  too  much  or  too  little  force  is  applied  to  the  underlying  soil 
and,  correspondingly,  either  too  little  or  too  much  force  is  applied  to  the 
structure. 

One  strategy  to  minimize  the  errors  from  the  Clough  and  Duncan  procedure 
is  to  use  as  few  excavation  stages  as  possible  in  the  simulation.  However,  this 
strategy  is  not  appropriate  for  problems  with  significant  nonlinearity  or  for 
excavation  problems  in  which  struts  or  tiebacks  are  placed  as  excavation 
proceeds. 

Both  the  method  based  on  the  accumulation  of  forces  and  the  method  based 
on  hybrid  finite  elements  satisfy  force  equilibrium  and  uniqueness  of  solution  for 
excavation  in  elastic  media.  However,  both  methods  are  complicated  and  difficult 
to  implement  in  a  finite  element  program,  and  the  requirements  of  force  equi¬ 
librium  and  uniqueness  of  solution  are  satisfied  only  at  the  expense  of  additional 
computational  effort. 

Of  the  four  general  methods  reviewed,  the  most  appropriate  method  for 
implementation  in  SOILSTRUCT-ALPHA  is  the  method  of  force  residuals.  This 
procedure  is  based  on  the  balancing  of  the  residual  forces  from  the  differences 
between  the  external  loads  and  the  equivalent  internal  loads  in  a  loaded  body. 

The  method  was  shown  to  satisfy  uniqueness  of  solution  and  force  equilibrium  in 
example  calculations.  It  was  also  noted  that  the  method  of  force  residuals  is  now 
the  most  widely  used  and  accepted  procedure  for  excavation  modeling.  Other 
advantages  of  the  method  of  force  residuals  are  that  it  uses  standard  finite  ele¬ 
ment  calculations  and  can  be  easily  extended  to  nonlinear  problems  in  a  con¬ 
sistent  manner.  Based  on  these  considerations,  it  was  recommended  to  update  the 
excavation  algorithm  in  SOILSTRUCT-ALPHA  by  using  the  method  of  force 
residuals. 

A  detailed  discussion  of  the  issues  related  to  implementation  of  the  method 
of  force  residuals  in  SOILSTRUCT-ALPHA  was  provided.  These  issues  include 
1)  the  type  of  finite  element(s)  to  be  used,  2)  the  implementation  of  interface 
elements,  3)  the  use  of  the  APHA  method  for  nonlinear  problems,  and  4)  issues 
related  to  the  computer  programming  of  the  method  of  force  residuals  in 
SOILSTRUCT-ALPHA. 

The  following  specific  recommendations  were  given  concerning  implemen¬ 
tation  of  the  method  of  force  residuals  in  SOILSTRUCT-ALPHA: 

a.  With  regard  to  the  finite  element  formulation,  it  is  recommended  to  keep 
the  current  QM5  finite  element  with  constant  stresses  in  the  code.  This 
has  the  advantage  of  retaining  the  improved  bending  performance  of  the 
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QM5.  It  is  also  suggested  to  implement  a  new  four-noded  isoparametric 
Q4  finite  element,  which  permits  the  calculation  of  element  stresses  at 
four  Gaussian  quadrature  points.  An  isoparametric  element  will  give 
better  resolution  of  stresses  and  excavation  forces  at  sharp  excavation 
boundaries  than  the  QM5.  Future  extension  to  SOILSTRUCT- ALPHA 
should  also  consider  higher  order  elements  such  as  the  eight-noded 
isoparametric  Q8  finite  element. 

b.  It  was  shown  that  interface  elements  can  be  accounted  for  in  similar 
manner  as  continuum  elements  in  the  method  of  force  residuals.  Example 
problems  were  used  to  show  the  significance  of  interface  behavior  in 
excavation  modeling. 

c.  It  was  shown  that  the  ALPHA  method  for  nonlinear  materials  can  be 
used  in  conjunction  with  the  method  of  force  residuals.  It  is 
recommended  that  excavation  be  performed  in  stages,  taking  out  a  few 
elements  at  a  time  (via  a  layer-by-layer  excavation  for  instance)  and 
applying  the  ALPHA-method  for  each  excavation  stage.  This  method  is 
preferred  because  both  geometric  and  material  nonlinearity  can  be  better 
accommodated. 

d.  Two  algorithms  for  implementing  the  method  of  force  residuals  in 
SOILSTRUCT- ALPHA  are  described.  Flowcharts  are  given  for  these 
two  algorithms. 
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Appendix  A 

Internal  and  External  Nodal 
Forces  Using  Isoparametric 
Elements 


To  clarify  the  nature  of  the  equivalent  external  nodal  forces  from  surface 
tractions  {Fexl }  (shown  in  Figure  2.3  and  given  in  Equation  3.3)  and  the  equiva¬ 
lent  internal  element  forces  {Fmt}  from  the  element  internal  stresses  (given  in 

Equation  3.4),  isoparametric  elements  will  be  used.  Isoparametric  elements  use 
the  same  shape  functions  to  describe  the  geometiy  (or  shape)  of  the  element  and 
the  unknowns  that  are  being  solved.  Further  details  of  the  isoparametric  formu¬ 
lation  can  be  found  in  finite  element  textbooks  (e.g.,  Cook,  Malkus  and  Plesha 
1989). 


A.1  External  Nodal  Forces  from  Surface  Tractions 

Consider  the  line  element  or  the  edge  of  a  two-dimensional  (2-D)  element 
•  shown  in  Figure  A.  1 ,  which  has  linear  displacements  u  and  under  the  surface 

tractions  (Tv.  To  represent  the  displacements  u  along  the  line  or  edge  of  an  ele¬ 
ment,  a  linear  function  of  the  coordinate  x  can  be  assumed: 

m(x)  =  a,  +  a2x  (A.l) 

The  values  of  the  interpolation  constants  a,  and  a2  can  be  obtained  from  the 
known  values  of  u  at  the  element  endpoints  at  x  =  Xj  and  x  =  x2 : 


Wj  —  +  a2x] 

u2-ax+  a2x2 


(A.2) 


Solving  these  two  equations  simultaneously  gives  the  values  of  the  two 
unknowns,  ax  and  a2 : 
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Figure  A.l.  One-dimensional  isoparametric  element  showing  (a)  local  coordinates  and  (b)  linear  function 
approximation 


«i~«t 


M,  -U. 

a2=~ - 1 

x,  “X 


Substituting  Equation  A.3  in  Equation  A.l  yields 
u(x)  =  N,(x)u,  +  N2(x)u2  =  [iVO)]{w} 

[^w]=[iv,w  tf2(*)],  {«}={"' j 

Af,(x)  =  l-^L,  N2(x)  =  ^±- 

X>  -X,  X,  -  X, 


(A.3) 


(A.4) 


where  the  Nl  values  are  called  element  shape  functions.  The  shape  functions  AT 
can  be  also  expressed  in  terms  of  the  local  coordinate  (5): 

u(s)  =  [N(s)]{u} 


n2(s)=1-^- 


5  =  2 


x-x, 


\X2-x]  J 


(A.  5) 


The  local  coordinate  5  takes  a  value  between  -1  at  x  =  xi  and  +1  at  x  =  x2,  and 
5  =  0  at  the  midpoint  between  x\  and  x2  (Figure  A.l). 


A2 
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Equation  3.3,  which  gives  the  equivalent  nodal  forces  { Fext ]  from  the  surface 
traction  a5,  is  actually  calculated  from  the  work  WT  done  by  the  surface  traction  a 
in  producing  the  surface  deformations  u: 


wt  =  \X*  Fex,(x)u(x)dx=  Fext  (5)w(5)c/5 


(A.6) 


Using  the  linear  distribution  of  displacement  u  in  Equation  A.4a, 


wt = r,1  f~'(s) 


1-5  1+5 


>ds 


(A.7) 


In  the  case  of  the  constant  distribution  of  traction  force  <f  along  the  line  element, 
as  shown  in  Figure  2.3a, 


Fext(s)  =  <50 


(A.8) 


f +1 


1-5  1+5 


\ds 


J 


(A.9) 


Performing  the  integration. 


>i{:> 


(A.  10) 


results  in  the  external  force  vectors  shown  in  Figure  2.3a 


fe}=^[ i  i] 


(A.  11) 


where  L-xt-xx  is  the  length  of  the  element. 

In  the  case  where  the  traction  force  CTV  varies  linearly  along  the  side  of  the 
element,  as  shown  in  Figure  2.3b,  the  same  displacement  function  [A(5)]  can  be 
used  to  express  the  variation  of  the  surface  traction  <f,  as 


^,(s)  =  [A,(5)  N2(s)]\J 


(A.  12) 


where  o,  and  are  the  magnitudes  of  the  surface  tractions  at  nodes  /  and  J. 
Substituting  in  Equation  A.7: 
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Performing  the  integration  yields  the  external  load  vector  shown  in  Figure  2.3b. 


(A.  14) 


Similar  procedures  can  be  derived  for  the  external  forces  from  body  loads  as 
given  in  Equation  3.3. 


A.2  Equivalent  Internal  Nodal  Forces  from 
Element  Stresses 

The  meaning  of  the  internal  force  vector  {F^,}  given  in  Equations  3.2  and 
3.5  will  be  illustrated  in  this  section  using  the  four-noded  isoparametric  quadri¬ 
lateral  element.  This  element  is  shown  in  Figure  A.2  in  the  global  {x,y}  coordi¬ 
nates.  Using  the  isoparametric  formulation,  the  element  can  be  mapped  to  the 
local  (or  natural)  {5,/}  coordinate  system  shown  in  Figure  A.2. 


Figure  A.2.  Quadrilateral  finite  element  (a)  in  the  global  coordinate  system  and  (b)  in  the  local  (natural) 
coordinate  system 
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Using  shape  functions,  the  horizontal  and  vertical  displacements  u  and  v  are 
interpolated  from  the  x  and  y  displacements  {«}  =  {w,  u2  u3  u4 }  and 

{v}  =  {vj  v2  v3  v4J  at  the  four  nodes  of  the  elements 


u  =  iVjW,  +  N2u2  +  N3u3  +  N4u4 
v  =  Nyx  +  N2v2  +  N3v3  +  N4v4 


(A.  15) 


In  the  isoparametric  formulation,  the  same  functions  used  to  express  the 
unknowns  u  and  v  are  used  to  express  the  shape  of  the  elements: 


x  =  N,  jc,  +  N2x2  +  N3x3  +  N4x4 


y  =  +  N2y2  +  N3y3  +  N4y4 


(A.  16) 


where  {jc}  =  {x,  x2  jc3  x4]  and  {7}  ={y,  y2  y3  y4}  are  the  coordinates 
of  the  four  nodes. 

Using  a  linear  combination  of  the  1-D  shape  function,  Ni ,  given  in 
Equation  A.5  for  both  the  local  axes  s  and  t,  the  four  2-D  shape  functions 
required  in  Equations  A.  15  and  A.  16  can  be  written  as 

,,  _(i-*Xl-0  Ar_(  1+sXl-O 

4  2~  4 

(A.  17) 

„_(l  +  sXl  +  0  Ar_(  l-s)(l  +  t) 

3“  4  4_  4 


The  [B]  -matrix  needed  in  Equation  3.3  to  determine  the  internal  nodal  forces  is 
defined  as 


[5] 


dNj_ 

dx 

0 


dN, 

dy 


dN,  d N, 
dy  dx 


(A.  18) 


The  derivatives  dN,  /dx  and  dN,  /dy  in  the  above  equation  can  be  obtained  from 
the  derivative  dN,  /ds  and  dNJdt  via  the  following  transformation: 
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where  [7]  is  the  Jacobian  matrix: 
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(A.20) 


From  Equation  A.  17,  the  derivatives  dN,  Ids  and  dN,  /dt  can  be  obtained  simply 
as 


dN. 

(1-0 
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4 

dt 

4 

li 

(l+o 

dN< 

_  0+0 

ds 
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In  the  local  5  and  t  coordinate  system,  the  equivalent  internal  force  vector  is 
expressed  as 


K,}=  JWM  dV  =  \+\  0*]r{o}|./|dHft 

V 

where  |j|  is  the  determinant  of  the  Jacobian  matrix: 


(A.22) 


dx  dy 
ds  ds 
dx  dy 
dt  dt 


(A.23) 
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The  derivatives  dx/ds,  dx/dt,  dy/ds  and  dy/dt  can  be  obtained  from 
Equation  A.  15  as 
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dx  _  ■A  dNj  dx  _  X-  dNt 
ds  ds  dt  dt 

dy  dN,  ,  dy  _y  dN, 

ds  tT  ds  y'  dt  £  dt  y 


(A.24) 


For  illustration,  only  the  value  of  the  horizontal  force  at  node  1  will  be 
calculated.  From  Equation  A.23, 


(A.25) 


Expanding  the  above  equations,  it  can  be  shown  that  Bu  and  |  J\  take  the 
following  values: 


Bn  -t2)+^(t3 -y,)+t{y2  -t3)] 


(A.26) 


where 
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1  -t 

t-s 

5-1 

t- 1 
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-(5  +  0 

s-t 

-(1  +  5) 
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1  +t 
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(A.27) 


(A.28) 


For  Element  2  in  the  example  problem  solved  in  Chapters  2  and  3, 

{*,}  =  {1  2  2  1}  and  {+,}={!  1  2  2}  .  Also,  in  the  same  example 

problem  in  Chapters  2  and  3,  the  stresses  are  constant  in  the  four-noded  isopara¬ 
metric  element.  Thus,  the  shape  functions  can  be  evaluated  anywhere  in  the 
element.  For  convenience,  the  local  coordinates  of  the  center  of  the  element, 
s  =  0  and  t  =  0  are  used.  These  values  give  |  J\ =1/4  and 


A,=i^[1+°(0)+°(-i)]=^ 


(A.29) 
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(A-30) 


Note  that,  in  general,  the  equivalent  internal  force  vector  given  in  Equa¬ 
tion  A.22  is  evaluated  by  numerical  integration.  The  most  common  technique 
used  is  the  Gaussian  quadrature.  In  this  technique,  the  function  to  be  integrated  is 
evaluated  at  some  optimum  locations  (called  quadrature  points),  the  value  of  the 
function  multiplied  by  a  weighting  factor,  and  summed  for  all  quadrature  points. 
In  similar  manner,  the  element  stiffness  [A:]  given  in  Equation  A.31  below  is  also 
integrated  numerically. 

(A.31) 


For  a  four-noded  isoparametric  element,  exact  numerical  integration  is 
achieved  with  four  quadrature  points,  while  for  an  eight-noded  element,  eight 
quadrature  points  are  required.  Consequently,  values  of  equivalent  internal  nodal 
forces  {/?„,}  and  the  element  stiffness  matrix  [*]  depend  upon  the  characteristics 

of  the  element  shape  functions  being  used  and  how  the  numerical  integration  is 
carried  out  in  Equations  A.22  and  A.31.  Thus,  it  is  not  sufficient  to  calculate 
excavation  forces  from  equivalent  internal  forces  alone.  To  ensure  consistency  in 
the  method  of  force  residuals,  the  following  conditions  should  also  be  met: 

a.  The  same  shape  function(s)  used  to  formulate  the  [5]  -matrix  for  the  load 

vector  in  Equation  A.22  should  be  used  in  the  [5]  -matrix  for  the  stiffness 
matrix  in  Equation  A.3 1. 

b.  The  equivalent  internal  nodal  forces  should  be  evaluated  using  the 
stresses  at  the  same  Gaussian  points  used  to  evaluate  the  stiffness  matrix. 


Condition  b  also  implies  that  the  same  level  of  Gaussian  quadrature  (i.e.,  the 
same  number  of  quadrature  points)  is  used  in  evaluating  the  stiffness  matrix  and 
the  equivalent  internal  load  vector. 
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Appendix  B 
Notation 


M 

at,  b, 

W 

[^] 

[*>] 

F, 

{FJ 


{4,} 


w 

[*] 

Kn,Ks 

Ka 

L 

w 

w 

[N] 


Interpolation  coefficient  matrix 
Interpolation  constants 

Self- weight  matrix  =  {0  J }T 

Strain-displacement  matrix 

Constitutive  matrix 

Nodal  force  at  node  i 

External  load  vector 

Internal  load  vector 

Element  stiffness  matrix 

Global  stiffness  matrix 

Interface  normal  and  shear  stiffnesses 

Coefficient  of  earth  pressure  at  rest 

Length  of  an  element  edge 

Coordinate  matrix 

Coordinate  matrix 
Element  shape  functions 
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[*] 

s,t 

U,  V 

M 

W 

Y 

{e} 

V 

M 


<*y 


Residual  force  vector 
Local  coordinates 

Horizontal  and  vertical  displacement 

Displacement  increment  vector 

Weight  of  an  element 

Unit  weight  of  soil 

Element  strains 

Poisson’s  ratio 

Element  stresses 

Surface  tractions 

Stresses  along  the  x-  and  y-axes 

Angle  of  friction 
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